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ABSTRACT 

We present an estimate of the black hole mass function (BHMF) of broad line quasars (BLQSOs) 
that self-consistently corrects for incompleteness and the statistical uncertainty in the mass estimates, 
based on a sample of 9886 quasars at 1 < z < 4.5 drawn from the Sloan Digital Sky Survey. We find 
evidence for 'cosmic downsizing' of black holes in BLQSOs, where the peak in their number density 
shifts to higher redshift with increasing black hole mass. The cosmic mass density for black holes 
seen as BLQSOs peaks at z ~ 2. We estimate the completeness of the SDSS as a function of black 
hole mass and Eddington ratio, and find that at z > 1 it is highly incomplete at Mbh ^ 1O 9 M0 and 
L/LEdd % 0.5. We estimate a lower limit on the lifetime of a single BLQSO phase to be tsL > 150±15 
Myr for black holes at z = 1 with a mass of Mbh — 10 9 AfQ, and we constrain the maximum mass 
of a black hole in a BLQSO to be ~ 3 x 1O 1O M0. Our estimated distribution of BLQSO Eddington 
ratios peaks at L/LEdd ~ 0.05 and has a dispersion of ~ 0.4 dex, implying that most BLQSOs 
are not radiating at or near the Eddington limit; however the location of the peak is subject to 
considerable uncertainty. The steep increase in number density of BLQSOs toward lower Eddington 
ratios is expected if the BLQSO accretion rate monotonically decays with time. Furthermore, our 
estimated lifetime and Eddington ratio distributions imply that the majority of the most massive 
black holes spend a significant amount of time growing in an earlier obscured phase, a conclusion 
which is independent of the unknown obscured fraction. These results are consistent with models 
for self- regulated black hole growth, at least for massive systems at z > 1, where the BLQSO phase 
occurs at the end of a fueling event when black hole feedback unbinds the accreting gas, halting the 
accretion flow. 
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1. INTRODUCTION 



Understanding how and when supermassive black holes (SMBHs) grow is currently of central importance in ex- 
tragalactic astronomy Observations have established correlations b etween SMBH mass and host g alaxy spheroidal 
properties, such as luminosity (e.g., Kormendv fe Richstqnel 119951; iMcLure fe Dunlopl 120011 120021), stellar velocity 
dispersion (Mbh~p relations hip, e.g., Gebhardt et al.l 20001: Merritt fe Ferraresell2001l : i Tremaine et al.ll2002D. concen- 
tration or Sersic .index Ce.g.. IGraham et al.l l2001t iGraham fe Driver! 12007), bulge mass (e.g., iMagorrian et al.llT99"p ; 



iMarconi fc Huntl 120031 : lHaring fe Rbd 120041 ). and binding energy (e.g. JAller fe Richst onc 20071 IHopkins et al.ll2007d) 
These correlations imply that the evolution of spheroidal galaxies and the growth of SMBHs are intricately tied to- 
gether, where black holes grow by accreting gas, possibly fueled by a major merge r of two gas-rich galaxies, until 
feedback energy from the SMBH expels gas and shuts off the accre tion process (e.g.. iSilk fc Reed 1998; Fa bianlH999t 
iBegelman fc Nathl 120051: iMurrav et all 120051: IHopkins et all l2009bl ) . This 'self-regulated' growth of black holes has 
been successfully appl ied in smoothed particle hydrodynamics simulations (|Di Matteo et al.ll2005l : iSpringel et al.ll2005l : 
Uohansson et aL | l2009D. and ha s motivated numerous models linki ng the SMBH growth, the quasar phase, and galaxy 
evolution (e.g.-lHaehnelt et al.lll998HKauffmann fc Haehneldl200tt lHaehnelt fc Kauffmannll2000l: iWvithe fc Loeblhooa 
Volonteri et al.ll2003t iCattaneo et al.ll2005t iDi Matteo et aljl2008t ISomerville et al.ll2008HHopkins et al.ll2008at iCrotonl 
20091 : ISijacki et aLl l2009: Bo oth fc Sch avc 2009; Shcn 2009, and references therein). Within this framework, the broad 



line qu asaiB phase persists after feedback energy from the black hole 'blows' the gas away (e.g.. IHopkins et al"]l2005aL 
2006b). The broad line quasar phase is e xpected to persist until the accretion rate drops low enough to switch to a 
radiatively inefficient accretion flow fe.g.. iChurazov et al.ll2005t IHopkins et al.ll2009al) . 

While major-mergers of gas-rich galaxies may fuel quasars at high redshift, and grow the most massive SMBHs, 
alternative fueling mechanisms are likely at lower redshift an d fainter luminositie s. Mergers alone do not appear to be 
sufficient to reproduce the number of X-ray faint A GN (e.g.. iMarulli et afll2007t ). and accretion of ambient gas (e.g., 
iCiotti fc Ostrikerll200ll : IHopkins fc Hernquistl 120061 ) . may fuel these fainter, lower Mbh AGN at lower z, resulting in 
an alternative growth mechanism for these SMBHs. Indeed, many AGN are observed to live in late-type galaxies out 
to z w 1 (e.g., iGuvon et al.1 120061 : iGabor et al.ll2009| ). and the X-ray luminosity function of AGN hosted by late-type 
galaxies suggests that fueling by mino r interactions or interna l instabilities represents a non-negligible contribution to 
the accretion history of the Universe ((Ge orgakakis et al. 2009). Furthermore, feedbac k from the SMBH may continue 
to affect its environment long after its growth through so-called 'radio mode' feedback (jCroton et al . 2006; B ower et al.l 
2006; Sijacki ct al. 2007). Observations qualitatitively support a model where the fueling mechanism for black hole 
growth depends on the mass of the host dark matter h alo, but regardless of the fueling mechanism , black hole feedback 
and accretion follow a similar evolutionary path (e.g.. iHickox et all 2009; C onstantin et al.ll2009D . 

Observationally, a signific ant amount of recent work has utilized the argument oflSoltanl (I1982D to indirectly map the 
growth of all SMBHs (e.g.. iSalucci et al.lll999l: lYu fc Tremaind 12001 iShankar et al.l 12004 12001 IMarconi et all l200l 
lYu fc Lull200l IHopkins et al.M2007bt IMerloni fc Heinzll2008D . Work along this line has used the correlations between 
Mbh and host galaxy spheroidal properties to infer the local distribution of Mbh for inactive black holes, which are 
assumed to be the relics of past AGN activity. The distribution of Mbh as a function of redshift is then estimated by 
stepping backward from the local distribution of Mbh, employing a conti nuity equation describing the 'flow' of black 
hole number density through bins in Mbh fe.g.. ISmall fc Bland ford 19921). The quasar luminosity function is used as 
a constraint on the rate of change in the SMBH mass density, because it traces the accretion of matter onto black 
holes, modulo the accretion efficiency and the bolometric correction. From this work, it has generally been inferred 
that black hole growth is dominated by periods of near Eddington accretion, with the most massive SMBHs growing 
first, and that many SMBHs have non- zero spi n. 

An alternative to the tech nique of ISoltanl (fl982l) for estim ating the SMBH mass function has been used by 
Siemiginowska 'fc Elvisl (|1997 | ) and lHatzi minaogl ou et al.l (|2001| ). These authors used a model for thermal- viscous 
accretion disk instabilities (jSiemiginowska et al.lll996l ) to calculate the expected distribution of luminosity at a given 
black hole mass. Based on this calculated distribu t ion, t hey use the quasar luminosity function to constrain the quasar 
black hole mass function. iSiemiginowska fc Elvisl (|1997| ) found evidence for black hole 'downsizing', with the peak of 
the quasar mass function shifting toward lower masses at lower redshift. 

Correlations between the SMBH mass, width of the broad emission lines, and lu minosity of the quasar contin- 
uum have made it poss ible to estimate Mbh for broad line quasars (BLQSOs) (e.g., iVestergaard fc Peterson! 120061 : 
iKellv fe~B cchtold 2007) , albeit with considerabl e statistical uncertainty of ~ 0.4 dex and various systematic effects 
(e.g.. lKrolikn200ll: IVestergaard fc Peterson! l200l iGreene fc Holl200l IMarconi et al.ll2008t iDennev et al.l [2009h . This 
offers an alternative to estimating SMBH mass functions based on the lSoltanl (|1982f l argument, because the mass func- 
tion may be estimated directly at all redshifts, and because the distribution of quasar emission line widths provides an 
additional observational constraint on the mass function. Estimates of Mbh obtained from the broad emission lines 
have been used to estimate the distribution of quasar black hole masses and Eddington ratios at a variety of redshifts 
(e.g.. IMcLure fc punlopll200l lVestergaaro1l2004l: iKollmeier et al.ll2006l: iNetzer fc Trak htenbrotil2007l: iShen et al.l l2008l: 
iFine et al.ll2008HKellv et al.ll2008HTrump et al.ir2009HLabita et al.ll2009af)~ 



The BLQSO black hole mass function (BHMF) maps the comoving number density and evolution of active super- 
massive black holes contained within broad line AGN, and is therefore a complete census of the population of these 
SMBHs over cosmic time. The BLQSO BHMF is important for a number of reasons, including the following: 

8 Throughout this work we will use the terms quasar and AGN to refer generically to broad line AGNs. No luminosity difference between 
the two is assumed, but they are assumed to have broad emission lines along our line of sight. 
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At hi gh redshift SMBHs with masses Mbh ~ 10 9 are already in place by z ~ 6 (e.g. jFan et al.ll200lUlJiang et al.l 
I2007D . and therefo re the high-z BLQSO BHMF places important constr aints on the formation and growth of 
seed SMBHs (e.g.. iVolonteri et al.ll2003l 12001 Eodato fc Nataralan1[200l . 

• If, after a fueling event, the growth of the SMBH persists until it becomes massive enough such that feedback 
energy begins to unbind the gas, the active SMBH will be seen as a BLQSO s hortly before entering: quies cence, 
and its fractional mass growth will not be significant during this time period ([Hopkins fc Hernqui st 2006) . The 
BLQSO BH MF thus (1) is re la ted to the distributio n of spheroidal binding energies in the central regions of the 
galaxy (|Hopkins et alJl2007al|cl : lYounger et al.ll2008l ). and (2) gives the nearly instantaneous increase in the AGN 
relic SMBH mass density. 

• Because mass is a fundamental physical quantity of SMBHs, we can use the BLQSO BHMF to estimate the duty 
cycle for broad line quasar activity as a function of mass by comparing the number density of all SMBHs with 
those seen as BLQSOs. The duty cycle can then be converted into an estimate of the lifetime of BLQSO activity, 
which is of significant importance for understanding the origin of BLQSO activity. This cannot be done using 
the quasar luminosity function. 

• The distribution of luminosi ties at a g iven BLQSO SMBH mass depend s on the quasar lightcurve (|Yu fe Lull2004t 
IHonkins fc HernauistJ[2006T : lYu fc Lull20M [Hopk ins fc Her nauist 2009) , which is a function of both evolution in 
the rate at which fuel is supplied to the accretion disk, a nd the time-dependent behavior of the accretion disk 
(jSiem iginowska & Elvis 1997t lHatziminaoglou et al.ll2001| ). Thus, understanding the distribution of L at a given 
MbHi or alternatively the distribution of Eddington ratio, gives insight into the BLQSO fueling mechanism and 
accretion physics. 

The BLQSO BHMF therefore provides an important observational constraint on models of SMBH growth, the onset 
and duration of quasar activity, quasar feedback, and galaxy evolution. 

There have been several estimates of the BHMF of BLQSOs ca l culated directly from the mass estimates derived from 
the broad emissio n lines ([Wang et al.ll2006t iGreene fc Hoi 120071: IVestergaard et all 120081: iVestergaard fc Osmer|[200l 
iKellv et aLll2009bl hereafter KVF09). In particular. IVestergaard fe Osmerl (|2009l ) found evidence for cosmic 'downsiz- 
ing' of black hole mass, in that the most m assive SMBHs are m o re common at high redshift, consistent with previous 
work on mapping black hole growth (e.g.. iMarconi et al.l [2004: Mc rlonil 120041: iShankar et alJl200l iMerloni fc Heinzl 



2008) , conclusions based on the quasar luminosity function (e.g.. iSteffen et al.l 120031: lUeda et al.l l2003t iCroom et al 



2004 lLa Franca et al.l 120051: i Hasinger et all 120051: iHopkins et all l2007bt iSilverman et al.l 120081 ) and the local active 
BHMF ( |Heckman et alf 2004) . However, a major concern with previous estimates of the BLQS O BHMF is uncorrected 
incompleteness a nd the additional broadening caused by the statistical errors in mass estimates (|Kellv fc Bcchtoldl l2007l : 
iShen et al.l[2008L KVF09). Because the massive end of the BLQSO BHMF falls off with increasing Mbh, the intrinsic 
uncertainty scatters more sources into higher Mbh bins than lower ones, potentially having a significant effect on 
the estimated number density of the most massive SMBHs. Furthermore, even if a sample is complete in luminos- 
ity, it is not necessarily complete in Mbh, and the completeness in Mbh depends on the unknown Eddington ratio 
distribution. Motivated by these issues, and the fact that the importance of the BLQSO BHMF demands that its 
determination be statistically rigorous, KVF09 developed a Bayesian statistical technique for estimating the BLQSO 
BHMF that self-consistently corrects for the incompleteness in Mbh and the statistical u ncertainties in the estima tes 
of Mbh- In this work we apply the technique of KVF09 to the SDSS BLQSO sample of IVestergaard et al.1 (|2008l ) in 
order to estimate the black hole mass function of SMBHs that reside in BLQSOs, and discuss the implications for 
SMBH growth, quasar lifetimes, and Eddington ratio distributions. 

2. THE DATA 

Our sample is drawn fr om t he Sloan Digital Sky Surv e y (SDSS) Data Releas e 3 (DR3) quasar sample as presented 
bv IRichards et al.l (|2006f ) and IVestergaard et al.l (|2008l ). IRichards et"afl (12006ft used 15,180 quasars from the SDSS 
DR3 to determine the quasar optical luminosity funct ion over 0.3 < z < 5. IVestergaard et aLl (|2008l) obtained black 
hole estimates for 14,434 of the quasars presented in IRichards et al.1 (|2006[ ) using scaling relationships between the 
width of the broad emis sion lines, continuum lumi nosity, and black hole ma ss. The details of the sample and fitting 
process are described in IRichards et all (120061) and IVestergaard et al.1 ([2008). For completeness, we briefly review the 
s pectral modeling used by IVestergaard et aLl <|2008f ) to extract the relevant quantities. 

IVestergaard et al.l (|2008l| modeled the observed quasar spectra using a power- law continuum, an optical-UV iron line 
spectrum based on I Zw I (IVestergaard fc Wilk es 2001: Veron-Cettv et al.ll2004l) . a Balmer continuum, and host galaxy 
template (jBruzual fc Charlotll2003D for objects at z < 0.5. The continuum luminosities used for the mass estimates 
are based on the power-law continuum fits. The emission lines used for the mass estimates in this work are Mg II and 
C IV, and were modeled using multiple Gaussian functions so to reproduce the line profile. Contributions from the 
narrow emission line region were subtracted from Mg II when appropriate, and line profiles with strong absorption 
were ignored. The end result of this analysis was a set of emission line FWHM and continuum luminosities, from 
which black hole mass estimates were calculated. 

We have performed a few additional cuts to the sample presented by IVestergaard et al.1 (|2008| ) before obtaining 
our final sample. First, we remove the sources located at z < 1. We do this because the quasar sample contains 
a significant number of extended sources below z < 1, and therefore their i-band magnitudes sometimes contain a 
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significant contribution from the host galaxy (see the discussion in Richards et al. (2006)). This reduces the effective 
flux limit at z < 1, creating an artificial second peak in the redshift distribution. While we could attempt to empirically 
correct the selection function to account for this, we find it easier to simply limit our analysis to 1 < z < 4.5. Finally, in 
order to make our analysis more robust against uncerta inty in the selection function, we omit any sources for which the 
value of the selection function of IRichards et all (|2006f) is less than 0.01, and force all values of the selection function 
to be zero that are < 0.01. After making these cuts, we were left with a sample of 9886 quasars. 

3. THE STATISTICAL MODEL AND DATA ANALYSIS 

We use an expanded version of the statistical model outlined in KVF09 to estimate the BLQSO BHMF. For complete- 
ness, we describe the important aspects of the technique developed by KVF09, and the reader is referred to KVF09 for 
further details regarding the technique and its derivation. Qualitatively, the technique of KVF09 assumes parameteric 
forms for the BHMF and the distribution of luminosities at a given Mbh- The average valu e of the broad line mass 
estima tes a t a given Mbh is fixed to be consistent with the scaling relationships reported by IVcstcrga ard fc Peterson! 
(2006) and iVestergaard fc Osmen <|2009f >: we assume that these scaling relationships produce unbiased estimates of 
Mbh- The BHMF and distribution of L at a given Mbh, hi combination with a selection function, imply an observed 
distribution of z, L, and FWHM. The technique of KVF09 attempts to recover the BHMF and distribution of L at a 
given Mbh by matching the observed distribution of z, L, and FWHM (or, equivalently, the mass estimates) that is 
implied by the model to the actual observed distribution. The posterior probability distribution is used to quantify 
how well the implied distributions match the observed distributions. As a result, the technique of KVF09 estimates 
the probability distribution of the BHMF, and of the distribution of L at a given Mbh, given the observed data set. 

3.1. Model for the Joint Distribution of Mbh, L, FWHM, and z 
We model the BLQSO BHMF as a mixture of K log-normal distributions: 



(j){\ogM B H,\ogz) 
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-1 K 



TT k 
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o(y-^fc) T E fe 1 (y-Aife) 
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where Ylk=\ ^ = ^ n ^ ms wor k we choose K = 4, based on our previous experience in working with simulated 
data sets (KVF09), and because we did not notice a significant difference in the results obtained using larger values 
of K. Here, N is the total number of BLQSOs in the Universe that could be seen by an observer on Earth at the 
time of the survey, y = (log Mbh, log z), is the 2-element mean vector for the A; th Gaussian functions, is the 
2x2 covariance matrix for the fc th Gaussian function, and x T denotes the transpose of x. In addition, we denote 
7r = (tti, . . . , ttr:)) /•* = ■ • ■ i P-k), and E = (Ei,...,E#). The variance in log Mbh for Gaussian function k is 
cr^ k = En.fc, the variance in logz for Gaussian function k is a z k = Y^22,k, and the covariance between log Mb// and 
logz for Gaussian function k is o m z,k — Ei2,fc. The parameters for the mass function are A, it, /i, and E. 

The distribution of luminosity density at a given black hole mass and wavelength A is assumed to also follow a 
mixture of J log-normal distributions: 



p(logL x \M B H) = 



J 

£ 

3=1 



2 ™h 



: exp 



1 f log XL \ - a j - a m j(log M B h - 9) 

2 I o-u 



(2) 



This represents an extension of the model of KVF09, which only used J = 1 log-normal distribution. We made this 
extension to incorporate additional flexibility in p(L\\Mbh), ensuring that the luminosity distribution, and therefore 
the black hole mass incompleteness correction, is robust to the particular parameteric form. For each log-normal 
distribution, the parameters for the distribution of L\ at a given Mbh are jj, aoj,a m j, and aij. We used J = 3 
log-normal distributions, as we did not notice a signficant change in the estimated values of p{L\\Mbh) when using 
J > 3. We further assess the robustness of our assumed form for p{L\\Mbh) in § 13.31 and show that our results should 
not be signficantly altered if the true form of p(L\\Mbh) is a power-law, as might be expected from some physical 
models for BLQSO lightcurves. 

In our analysis we use the luminosity density at 1350A in order to minimize bias at the highest redshifts introduced 
from extrapolating the power-law continuum beyond the spectral window. At z > 1.8 the rest frame A = 1350 A 
falls within the observed spectral window for the SDSS sources used in this work, while a smaller redshift window is 
available when using the luminosity densities calculated at other wavelengths. Furthermore, the z ~ 1 distributions 
of bolometric luminosity did not exhibit any significant difference when using the luminosity density at 1350A, as 
compared to that calculated using the luminosity density at A > 1350A, suggesting that significant biases at lower z 
are not introduced by extrapolating the power-law continuum. 

We can connect the parameters in Equation @ to the distribution of Eddington ratio and bolometric correction. 
The monochromatic luminosity is related to the Eddington luminosity ratio T EdcE\ and bolometric correction C\ as 



XL\ 



L3xl0 3 8 r^M B H r 



[erg s ]. 



(3) 



C x M Q 

9 We will occasionally use F Edd to denote the Eddington ratio, instead of L/L Edd . We do this in certain instances where it is more 
notationally convenient to do so. 
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Therefore, Equation implies that for the j th log-normal distribution we are assuming that on average log(r Edd/C\) — 
aoj — 47.11 + (a m j — 1) log Mbh, with a Gaussian scatter about this mean value of standard deviation aij. We do 
not make any formal attempt to prohibit Equation ([2]) from allowing values of L/L B dd > 1, as this would require 
us to make an assumption about the bolometeric correction. However, as we will show in § 14.51 our estimate for 
p{L\\Mbh) im plies only a negligible fractio n of BLQSOs with L/LEdd > 1> assuming a constant bolometric correction 
of C1350 = 4.3 (jVestergaard k Osmer|[200l . 

The distribution of emission line widths at a given luminosity density and black hole mass is modeled as a log-normal 
distribution: 

p{\ogFWRM BL \L BL ,M BH ) = 

1 f 1 (log FWHMgx - Af L + 1/4 log L BL - 1/2 log M BH ) 2 

V 27T ( a BL + °FWHm) I- 2 a BL + a FWHM 

Here, FWHMbl is the line width for a particular broad emission line, Lbl is the luminosity density used as an estimate 
for the broad line region size for a particular broad emission line, ctfwhm is the measured uncertainty in FWHM due 
to measurement error, and /3q L and <Jbl are the parameters for Equation Q for a particular broad emission line. We 
do not make any attempt to correct for radiation pressure on the broad emission lin e clouds dMarconi "et"ap2 008l as 
its importance and effects are currently poorly understood (e.g., for a discussion see iVestergaard k Osmerll2009D . 

Equation (0} implies that on average FWHM oc Ml'x/L 1 ^, or equivalently M BH oc L^FWHM^. Therefore, 
we can use the results obtained for the broad e mission line scaling estimates o f M B h to fix We use the mass 
scaling relationship for C IV that is presented by IVestergaard k Peterson! ((2006:), and the relationship for Mg II that 

is presented by IVestergaard k Osmerl (|2009l ). As noted in KVF09, these scaling relationships imply /3^f sl1 = 10.61 and 
/3q IV = 11.33, and we fix /?o to these values. 

In addition, the dispersion in FW HM B l at a given L B l and Mbh can be related to the scatter in the mass estimates 
based on the scaling relationships. Vesterga ard k Peterso n (2006) find the statistical uncertainty in the broad line 

1/2 

mass estimates to be 0.36 dex for C IV. Therefore, because FWHM oc M BH , ctbl(CIV) « 0.18 dex. Likewise, 
the intrinsic uncertainty in the broad line mass estimate for Mg II is found to be ~ 0.4 dex (Vestergaard et al., in 
preparation), and therefore ctbl (Mgll) rs 0.2 dex. However, there have been indications from previous analysis of 
flux limited surveys, which probe higher redshifts and a narrower range in luminosity than that exhibited by the 
objec ts with reverberation mapping data, that the intrinsic scatter in the mass estimate s may be smaller than rs 0.4 
dex (|Kollmeier et all 120061: iShen et all 120081: iFine et all 120081: iSteinhardt k Elvisll2010aHbT) . This may be caused by 
correlated scatter in the mass estimates with luminosity (e.g., IShen k Kellvl I201QT) or redshift, or a dependence on 
o~bl with luminosity or redshift. Both of these possibilities would decrease the d ispersion in the scatt er in the mass 
estimates when only probing a smaller range in L, or higher redshifts. In addition, iMarconi et al.l (|2008l ) find a smaller 
scatter in the masses estimated using H/3 when one corrects the reverberation mapped masses for radiation pressure. 
Because of the possibility for a smaller scatter in the mass estimates, we perform our analysis with both ctbl fixed to 
« 0.2 dex, and with <jbl as a free parameter. 

3.2. The Posterior Distribution and Fitting Technique 

The technique for estimating the BHMF developed by KVF09 takes a Bayesian approach for performing statistical 
inference, meaning that it calculates the probability distribution of the mass function, given the observed data. All 
the information regarding the BHMF and the parameters for the statistical model is contained within the posterior 
probability distribution, which is defined as the probability distribution of the model, given the observed data. KVF09 
derived the posterior distribution for the statistical model described in the previous section. Denoting the model 
parameters for the shape of the BHMF as 9 = (ir, fi, S, 7, a 0l ct m , 01), the posterior distribution is 

n 

p(#|FWHM, L, z) oc P {9) [p{I = l|0)p J]p(FWHM i5 L x>i , Zi \9), (5) 

i=l 

where the number of data points is n, p{9) is the prior on 9, and p(I = 1\9) is the probability as a function of 
9 that an BLQSO makes it into the SDSS DR3 catalogue. We note that if the dispersion in the mass estimates, 
o-bl, is also treated as a free parameter, then 9 also contains ctbl- The joint distribution of FWHM,i, and z, 
p(logFWHMi,logL.\,i,logZi|#), is given by Equations (30)-(40) in KVF09, modified to use the mixture form for 
p{L x \M B h)- We use'the the Mg II line at 1 < z < 1.6, both the Mg II and the C IV line at 1.53 < z < 1.6, and the 
C IV line at z > 1.6. We do not use H/3 because we limit our analysis to z > 1. At z ~ 0.8 H/3 is shifted into the 
water vapor bands, which tends to decrease the FWHM accuracy; Mg II is similarly affected at higher redshifts. In 
addition, we see systematic changes in the Mg II FWHM distribution above a redshift of 1.6, suggesting the presence 
of biases, which need further investigation (M. Vestergaard et al., in preparation). The posterior distribution for the 
BHMF normalization, given 9 and n, is given by Equation (16) in KVF09. 

The inclusion probability as a function of 9 is calculated by averaging the SDSS DR3 selection function over the 
distribution of L\ and z (see Eq.(46)-(49) in KVF09). In order to simplify our analysis we ignore the lower limit of 
FWHM > 1000 km s" 1 on the emission line width for the SDSS DR3 sample. The distribution of FWHM for the 
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SDSS DR3 falls off before reaching FWHM = 1000 km s _1 , and it does not appear that imposing the lower limit 
results in a non-negligible fraction of the BLQSO population being missed. Therefore, any correction for the lower 
limit in FWHM is negligible, and we ignore it. In this case, the inclusion probability is 



P(I = 1\0) = - 





f 45 s(L x ,z) 




=oJz=i L x zhil0 
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J2ljN(logL x \l 


fe=l 





\\hj(z), Vi^j) 



N(logz\n Ztk , a 2 k ) dz dL x (6) 
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Here, Q. = 1622 deg 2 is the effective sky area of the SDSS DR3 sample ^Richards et al.l 120061) . s(L x ,z) is the SDSS 
selection function, N(x\p, a 2 ) denotes a Normal distribution with mean p, and variance c 2 , as a function of x, p m ^k and 
p z ,k are the mean values of log Mbh and logz for the k th Gaussian function, respectively, and p m z,k is the correlation 
between log Mbh and logz for the k th Gaussian function. Note that lk(z) and Vi : k define the mean and variance in 
log La at a given redshift for the k th Gaussian function. The integral in Equation ^ is over 1 < z < 4.5 because we 
have removed the sources at z < 1, and there are no useable broad emission lines at z > 4.5. 

Equation © is in terms of the selection function with respect to the luminosity density. As mentioned above, we 
use the luminosity density at 1350A in this work. However, IRichards et al.l (120061) repo r t thei r selection function in 
terms of the i-band magnitude. We can convert the selection function of IRichards et al] (j2006l ) to be in terms of the 
BLQSO power-law continuum L1350 through the equation 

/oo 
s(i,z)p(i\Li 350 ,z) di, (9) 
-00 

where p(i|Li35o, z) is the distribution of «-band magnitude at a given L1350 and z, and s(i, z) is the SDSS DR3 selection 
function in terms of i and z. We model the distribution of i magnitudes at a given L1350 in different redshift bins as 
a student's t distribution with mean that depends linearly on logLi35o: 

r[(K*) + i)/2] ( 1 ( i - Mz) - Bjjz) iog£ 135 o A 2 \ ~ Kz)+1)/2 

The student's t distribution converges to the normal distribution for v — > 00; for finite v the t-distribution is more 
heavy tailed than the normal distribution, and we have found it to better describe the distribution of i-magnitudcs 
at a given Z/1350 and z. The parameters Ai(z), Bi(z), &i(z), and v(z) are fit by maximizing their posterior probability 
distribution, given the observed set of i mag nitudes and £1350- The posterior distribution is given by inserting the 
assumed distributions into Equation (40) oflKellvl (l2007f). and fo r simplicity we assume a simple flux limit of i — 19.1 
for z < 2.7 and i = 20.2 for z > 2.7 (e.g.. IRichards et al.|[2006D . At z < 2.7 the redshift bins used in Equation ([ID]) 
have width Az = 0.1, while at z > 2.7 the redshift bins have width Az = 0.3. 

In this work we constrain the parameters of our statistical model to be within certain limits, but in general assume 
a uniform prior on 8. Our prior is uniform on the parameters within the limits 44.5 < ao < 46.5, —1 < a m < 3, 0.1 < 
<ti < 2,7 < fi m _k < 10, log 1 < [i z _k < log 4. 5, 0.1 < <T m ^,o'z.k < 2, and —0.98 < p mz ,k < 0.98. We also assume 
a uniform prior on 7r and 7, subject to the elements of n and 7 summing to unity. The limits on a m were chosen 
because we did not consider it realistic that £1350 would depend on Mbh outside of the range of dependencies spanning 
£1350 oc I /Mbh to L1350 cx M BH , and the limits on ao were chosen to restrict the average value of L/LEdd to be 
within 0.01 < L/LEdd < 1 for BLQSOs with Mbh = 10 9 M Q , assuming a bolometeric correction of C1350 = 4.3. The 
limits on <j\ were chosen because A log L w 0.1 is comparable to the grid spacing on which the selection function was 
computed by [Richards et al.l (|2006[ ). and we did not consider it likely that the dispersion in £1350 at a given Mbh 
would be greater than 2 dex. The limits on the BHMF parameters were chosen so as not to extrapolate the BHMF very 
far beyond the detectable range of Mbh- As such, in this work we limit our analysis of the BHMF to Mbh ^ 10 7 Af Q . 

As mentioned before, there exists the possibility that, in the range of L and z we probe, the uncertainty in the mass 
estimates may be smaller than the commonly quoted ~ 0.4 dex. If the error in the mass estimates is correlated with 
luminosity, then the variance in the mass estimates at a given luminosity and mass is reduced to 

Var{\ogM BL \L) - Var(logM BL )(l - p% L ). (11) 

Here, Mbl denotes the broad line mass estimate, Var (log Mbl) is the variance in the mass estimates about the true 
mass, typically thought to be ~ 0.4 2 dex 2 , and pbl is the correlation with luminosity in the scatter in the mass 
estimates about the true mass. Because we probe a somewhat narrow range in L\, when we estimate the parameter 
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<7bl in Equation (|4|, we are really estimatin g 2a bl ~ y V ar (log Mb l\L), and we therefore need to also construct a 

prior for obl- We do this by first noting that iVestergaard fc Peterso n (2006) estimated the dispersion in the scatter in 
the broad line mass estimates about the reverberation mapping estimates using 27 AGN. Therefore, the appropriate 
prior distribution for Var(\og Mbl) is the posterior probability distribution of the variance in the mass estimates, 
given the reverberation mapping sample. Since we assume that the scatter in the mass estimates about the true mass 
is log-normal, the probability distribution of their variance follows using a standard result fro m Bayesian statistics, 
and is a scaled inverse \ 2 distribution with 26 degrees of freedom (e.g.. iGelman et all |2004 ). However, there are 
currently no constraints on the value of Pbl-, and we use a uniform prior on its value. Our prior on ubl is then 
calculated by combining these two priors according to Equation This results in a broad prior which peaks at 

V ar (log M B l\L) « 0.4 2 dex 2 , and falls off slowly to zero as V ar (log Mb — > and Var(\og Mbl\L) — > 0.6 2 . 

In addition to the above constraints, we also impose the prior constraint that the number density of BLQSO SMBHs 
must never exceed the local number density of all SMBHs. In principle, this constraint could be violated if a large 
number of 'wandering' black holes are present. These wandering black holes would be ejected from g alactic nuclei in 
the late stages of a merger due to asymmetric emission of graviational radiation (e.g.. IVolonte ri 2007). However, this 
constraint would only be violated if the binary black hole system is ejected after or during the broad line quasar phase. 
While this would certainly be a very intriguing result, we do not test it here; indeed, our results are unaffected by 
this prior constraint as the estimated BHMF is always below the local value for all random draws from the posterior 
probability distribution. 

As described in KVF09, we do not work with Equation ([5]) directly, but instead use a Markov Chain Monte Carlo 
(MCMC) sampler algorithm to obtain a set of random draws of 9, distributed according to Equation ([5]). Because the 
posterior distribution described by Equation ([5]) is multimodal due to ambiguities in labeling the Gaussian functions, we 
label the BHMF Gaussian functions in order of increasing implied mean flux, and the p(L\Mbh) Gaussian functions in 
order of increasing mean luminosity at Mbh = 10 9 A/©. In addition, in order to make our MCM C sampling algorithm 
robust against additional possible multimodality, we include parallel tempering (e.g.. lLivj|l2004fl in our Markov Chain 
Monte Carlo (MCMC) algorithm to facilitate sampling from the different modes. For each value of 9 we obtained 
from our MCMC random number generator, we also obtain a value of the BHMF normalization, N, by drawing from 
a negative binomial distribution with parameters n and p(I = 1\9) (KVF09). The random realizations of 9 and N 
then define a sample of random realizations of the BHMF obtained from the posterior probability distribution of the 
BHMF, given the observed set of mass estimates, luminosity densities, and redshifts. We ran our MCMC sampler for 
2 x I0 5 iterations each, keeping every 20 th iteration. 

3.3. Robustness to Incorrect Assumptions Regarding the BLQSO BHMF and Eddington Ratio Distribution 

Although we have assumed a flexible parameteric form for the BLQSO BHMF and p(L\Mbh), this form is incon- 
sistent with more physically- motivated BLQSO p(L\Mbh), such as migh t be expected from a power-law decay in the 
BLQSO accretion rate (e.g., lYu et all 120051 : Hopkins & Hernquist 2006], also see discussion in § 15.11) . Instead, our 
assumed log-normal mixture form for p(L\Mbh) was motivated by mathematical convenience, as some of the integrals 
necessary for evaluation of the posterior distribution can be done analytically (KVF09). If we were to use a form which 
required numerical integration, we would have to perform over ten billion numerical integrals in our MCMC sampler, 
which is computationally prohibitive. Moreover, we also used the mixture form to allow flexibility in the estimated 
p(L\Mbh), so that our assumed form should be able to approximate many different forms for p(L\Mbh)- 

In order to assess the impact of our chosen parameteric form on the inferred BHMF, we simulated a data set where 

the distribution of Edd ington ratios was assumed to have a power-law form p(T Edd\MBH) oc ^Edd' 1 ^^ with /3 = 2 
(see discussion in § I5.I|) . The distribution of bolometric corrections was log- normal with geometric mean C1350 = 5 
and dispersion of 0.2 dex, and the BHMF normalization was N = 2 x 10 6 ; all other aspects of the simulation were 
done in the same manner as described in § 6.1 of KVF09. 

The results are illustrated in Figure [TJ where we compare the true and estimated BHMF for the simulated sample at 
2 = 2, and the true Eddington ratio distribution with that inferred assuming the statistical model described in § 13.11 
Here, and elsewhere in this paper, in addition to a single 'best-fit' mass function, we also plot 100 random draws of 
the mass function from its probability distribution, generated by our MCMC random number generator. The spread 
and density of the random draws of the BHMF, and any quantities derived from it, give a visual representation of the 
uncertainty in these quantities, with the spread on the random draws constraining the BHMF. Because we use 100 
random draws, the probability of a quantity falling within a certain area on a plot can be estimated by counting the 
number of random draws that intersect that area. For this example, the BLQSO BHMF inferred assuming a mixture of 
J = 3 log-normal distributions for p(L\\Mbh) is able to recover the true BHMF and Eddington ratio distributions, at 
least when the true distribution of L/L<Edd is p{J^Edd) °c ^~Edd- Although this test is far from exhaustive, we consider 
it reasonable to conclude that our flexible form for the BHMF and p(L\\Mbh) is able to accurately approximate the 
true forms, and therefore our results are robust against errors resulting from using an incorrect parameteric form for 
these distributions. 

4. RESULTS 

4.1. Evaluating the Fit: How Uncertain are the Mass Estimates? 
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Fig. 1. — (a) True BHMF for the simulated sample described in § 13.31 having a power-law distribution of Eddington ratios over 0.01 < 
L/LEdd < 1 (thick red line), compared with the BHMF estimated assuming a mixture of J = 3 log- normal distributions of Eddington 
ratios (thin black lines); each of the thin s olid lines denotes a random draw from the probability distribution of the BHMF, given the 
observed data and assumptions outlined in § 13.11 Here, and in all figures in this work, we plot 100 random draws of the function of interest, 
so the probability that the BHMF, say, has a certain value in a given range can be estimated by counting the number of random draws of 
the BHMF that fall within that range. For this example, the BHMF estimated assuming a mixture of log- normal distributions for L/L^^ 
is able to recover the true BHMF, implying that our mixture form is robust against mispccification of the parameteric model, (b) True 
Eddington ratio distribution for the same simulated sample (thick red line) , compared with the estimated distribution assuming a mixture 
of J = 3 log-normal distributions (thin black lines). In this plot we have convolved the power-law form of the distribution of L/L^dd with 
the scatter in the bolometeric correction used in this simulation, to incorporate the error in the estimated Eddington ration distribution 
introduced from assuming a constant bolometeric correction. The mixture of log-normals form is able to adequately approximate the true 
distribution of L/L^dd- 
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Fig. 2. — Distribution of 2, £1350, and the mass estimates for our sample (histograms), compared with the distributions implied by our 
best-fit BHMF and distribution of L1350 at a given Mbh (red squares). In the top row of plots we fix the statistical error in the mass 
estimates to be 0.4 dex, while in the bottom row we allow it to be a free parameter. The error bars denote the la uncertainties on the 
implied distributions, although often they are smaller than the red squares. The statistical model described in § 13.11 is able to reproduce 
the observed distributions only if we allow the standard deviation in the statistical error in the mass estimates to be a free parameter, 
implying that the standard value of ~ 0.4 dex is too large for this luminosity and redshift range. 



We used our Bayesian method to derive the BLQSO BHMF from the SDSS DR3 quasar sample, both holding 
the magnitude of the scatter in mass estimates fixed to 0.4 dex, and treating the amplitude of the scatter as a 
free parameter. However, before discussing the results, we first evaluate how well the statistical model described in 
§ 13.11 fits our sample. In order to do this, we compare the distributions of redshift, luminosity, and broad line mass 
estimates of our sample to the distributions implied by our model, as described in KVF09. Figure [2] compares the 
observed distributions of z, \L\(1350A), and Mbl to those implied by the statistical model described in § 13.11 The 
implied distributions were calculated by first simulating a sample of Mbh and z from a BHMF randomly output from 
the MCMC sampler. Then, for each value of Mbh, we simulated a value of ALa(1350A) and mass estimate Mbl- 
Finally, we applied the selection function given by Equation ([9]) to the simulated data set. This was repeated for each 
realization of the BHMF obtained from the MCMC output, in order to account for the uncertainty in our estimated 
model parameters. 

We are not able to fit the mass estimate distributions if we keep the statistical scatter in the mass estimates fixed 
to 0.4 dex. In particular, this predicts a distribution of mass estimates that is too broad compared to the actual 
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distribution. In fact, the observed dispersion in the mass estimates for our SDSS sample is sa 0.35 dex, smaller than 
that expected even if all objects had the same mass. However, if we allow the dispersion in the mass estimate error 
to be a free parameter, we are able to obtain a good match to the data. Our best-fit value for the scatter in the mass 
estimates at a given mass is 0.18 ± 0.01 dex for Mg II and 0.13 ± 0.01 for C IV. Our result that the scatter in the 
mass estimates for high L and z must be smaller than ~ 0.4 d ex is consistent with what has been found in previous 
work (Ko llmeier et al.l [2006: Sh en et al.l l2008: Fi ne et al~l 12008). and our best-fit va lues of the standard deviati ons in 
the mass estimate errors are consistent with the upper limits recently calculated by iSteinhardt fe Elvisl (|2010bl ). 

The origin of this smal ler scatter i s uncl ear, a nd there may be sev eral different possibilities. One possibility, as 
mentioned earlier and by iShen et ail (|2008[ ) and iShen fc Kellvi ()2010[ ). is that the error in the mass estimates may 
be correlated with luminosity. If this is true, then an error of ~ 0.4 dex represents the error in the mass estimates 
when averaging over a broad range in L, as was done for the reverberation mapping sample, while a smaller scatter 
of ~ 0.15 dex represents the error in the mass estimates when one is limited to a more narrow range in L, as the 
SDSS quasar sample is. It is unclear why the error in t he mass estimates wou ld be correlated with luminosity, but one 
possible source unaccounted for is radiation pressure. iMarconi et al.l (|2008l ) argue that virial mass estimates should 
be corrected for radiation pressure. They find a correction that implies a steeper dependence on luminosity than 
Mbl oc L ' 5 , especially for sources with high L/LEdd- Under their model, if one does not correct for radiation pressure 
then one will tend to underestimate the mass with increasing L, producing a correlation between the error in the mass 
estimates with luminosity, and possibly producing a smaller s catter in the mass estimate errors over a narrow range in 
luminosity. It is interesting to note that IMarconi et all (|2008l ) find a smaller scatter in the mass estimates of ~ 0.2 dex 
when correcting for radiation pressure for the reverberation mapped sample, which covers a larger range in luminosity. 
However, more work is need to understand the importance of radiation pressure, and if it can produce the observed 
smaller scatter in the mass estimates over the range in luminosity we probe. 

Another possibility is that the error in the mass estimates may not be correlated with L, but the dispersion in the 
errors may decrease with increasing L or z. Unfortunately, the number of AGN with Mbh estimated from reverberation 
mapping is too small to test this, and dominated by sources at lower L and z. The third possibility is that the virial 
mass estimates are biased at high L and z, at least for Mg II and C IV, and are only marginally related to the actual 
masses. The R-L relationship is well established for H/3, and does not require a large extrapolation to the luminosities 
probed in our sam ple, and thus we do not expect H/3-based mass estimates to be significantly biased in this range 
(Vcstergaard 2009) . The situation is less clear for Mg II and C IV. There is only one reliable time lag for the Mg II 
emission line ()Metzroth et al.ll2006D. An R -L relationship has been estimated for the C IV line over a broad range in 
luminosity and redshift (jKaspi et al.ll2007|) . including those probed in our study, and the R-L relationship is similar 
for both H/3 and C IV. Unfortunately, the C IV R-L relationship is estimated from only eight data points, and further 
work is needed in order to understand the Mg II- and C IV-based mass estimates and their errors. 

Throughout the rest of this work will we focus on the results obtained from allowing the dispersion in the mass 
estimate error to be a free parameter, a value of ~ 0.4 dex is clearly ruled out. However, we note that we have 
performed the same analysis for both cases, and while the quantitative details change, the scientific conclusions are 
unaffected by treating the amplitude of the scatter as a free parameter. The only exception is that we infer a much 
more narrow distribution of Eddington ratio if we fix the amplitude of the scatter in the mass estimates to be ~ 0.4 
dex. In addition, the results reported in this section highlight the need for more reverberation mapping studies in 
order to better understand the nature of the errors in the mass estimates. 



4.2. The Black Hole Mass Function for Broad Line AGN 

Figure [3] shows the BLQSO BHMF at several redshifts. Our es timated BHMF is compa red with an estimate of 
the local mass function of all SMBHs, and t he BHMF reported bv Vestergaard et al.l (|2008l ). obtained from binning 
up the broad line mass estimates. Following Me rloni fe Heinzi (12008 ). the local BHMF was computed to be near the 
middle of the uncertainty range reported bv lShankar et all ((2009) by convolving a Schechter function with a log- normal 
distribution with standard deviation 0.3 dex, chosen to be cons istent with the scatter a bout the Mbh~o~ relationship; 
the parameters for the Schechter function are those reported by Me rloni fc Heind (|2008| ). Also, we note that early type 
galaxies dominate the local BHMF at M BH > 4 x 1O 7 M (|Yu fc Lull20~0ll . 

In order to focus on the region of the BHMF that is robust against uncertainties in the selection function, as well as 
against uncertainty on the Eddington ratio distribution, we estimate the black hole mass completeness for our sample 
as a function of z. Our best estimate of the SDSS completeness as a function of Mbh and z is shown in Figure IH 
and the 10% completeness limit is marked by a vertical line in Figure |3] The black hole mass completeness depends 
on both the completeness in luminosity, and the assumed distribution of L at a given Mbh- As can be seen from 
Figure SI at z > 2 the SDSS quasar sample is only rs 10% complete at Mbh ~ 10 9 M©, becoming more incomplete at 
lower masses. At masses much lower than Mbh ~ 1O 9 M , the estimated mass function almost completely depends on 
extrapolation from the set of BHMFs and Eddington ratio distributions that fit the observed data well, constrained 
by our assumed paramctcric forms. Therefore, we stress that below Mbh ~ 10 9 M Q the BLQSO BHMF must be 
interpreted with caution, and in this work we will try to focus on what we can infer from the high mass end of the 
mass function. 

We also estimate the completeness of our sample as a function of L/LEdd-, assuming a constant bolometric correction 
of C1350 = 4.3. The estimated completeness depends on the selection function and the estimated BHMF, as the selection 
function depends on luminosity, which is defined by the BHMF at a given L/LEdd- The Eddington ratio completeness 
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Fig. 3. — BLQSO BHMF (thin solid lines) obtained using our Bayesian approach, compared with the local BHMF fo all SMBHs (dashed 
line), and the BHMF from Vcstcreaard ct al. (2008, solid red line with points); as in Figure[T] each thin solid line denotes a random draw 
of the BHMF from its probability distribution. The thick green line is the median of the BHMF random draws, and may be considered 
our 'best-fit' estimate. The vertical line marks the mass at which the SDSS DR3 sample becomes 10% complete. 
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Fig. 4. — Estimated completeness in black hole mass (left) and Eddington ratio (right) for the SDSS DR3 quasar sample of Richards ct al. 
(2006), calculated using our assumed distribution of luminosity at a given black hole mass (Eq.[2]). The red, green, and blue lines denote 
the 10%, 50%, and 90% completeness levels, respectively. The SDSS sample is highly incomplete at Mbh $ 10 9 Mq and L/Lgdd S; 0.5. 
Note that the incompleteness at the highest masses is due to the upper flux limit of the SDSS. 



is also shown in Figure [4] The SDSS is only ~ 50% complete for sources radiating at the Eddington limit, and < 10% 
complete for sources radiating at < 10% of Eddington. This heavy incompleteness is due to the fact that, for our 
estimated BLQSO BHMF, half of BLQSOs have black holes that are not massive enough to make the SDSS flux limit 
even if they radiate at Eddington. We further discuss the Eddington ratio distribution in Section § 14.51 

The difference between the binned estimate of the BLQSO BHMF, calculated from the estimate of 
(2008), and our estimate, is the result of the difference in statistical methodology employed by 
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2008) and our work. First, the statistical erro r in the broad line mass estimates results in a broader inferred BHMF 



when one simply bins up these estimates (e.g.. lKellv fc Bec htold 2007, KVF09). Second, the 1/V a technique corrects 
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Fig. 5. — Evolution in the number density of BLQSO SMBHs for four different values of Mbh- A larger fraction of the most massive 
SMBHs are seen as BLQSOs at high er redshift, an effect com monly referred to as black hole downsizing. Downsizing of BLQSO SMBHs 
has also been seen in the LBQS (Vestergaard & Osmcr 2009). Symbols are as in Figure [3] The wiggles seen in the lowest mass bin are an 
artifact of fitting a Gaussian mixture model, unlikely to be real, and probably due to small unaccounted for errors in the selection function, 
which can become magnified due to the the fact that we are incomplete in this mass bin. 

for incompleteness in flux, but not black hole mass. As a result, the 1/V a corrections only partially correct for 
incompleteness in Mbh, an d the estimated binned BHMF will still suffer from incompleteness, especially at the 
low mass end. T he two effects combined result in a systematic shift in the estimated BHMF toward higher Mbh 
(|Shen et al.ll2008l KVF09). However, the Bayesian approach outlined in KVF09 is able to self-consistently correct 
for these two effects in a statistically rigorous manner, conditional on the survey selection function and assumptions 
implicit within t he statistical mode l. In s pite of these differences in methodology, our estimated BHMF agrees fairly 
well with that of [Vestergaard ct al. (200<|) at the high mass end, considering the differences in the methodology, and 
the two do not strongly diverge until values of Mbh where the SDSS is highly incomplete. The poorer agreement 
between the two estimates at z < 2 is due to the fact that the BHMF is estimated using the Mg II emission line in 
this redshift range, which we find to have a higher statistical error than that of C IV. As a result, the our correction 
to the BHMF due to the error in the mass estimates is greater at these redshifts. 

In Figure [5] we show the evolution in the comoving number density of SMBHs in BLQSO at four different masses. 
As is evident, the number density of higher mass BLQSO SMBHs peaks at higher redshift, where the number density 
of BLQSO SMBHs of M B h ~5x 10 8 M q peaks at z ~ 1.5, and the number density for M B h ~5x 10 9 M q peaks at 
z ~ 2.5. This result is consistent with what has commonly been referred to in the literature as black hole 'downsizing', 
w here the most massive SMBHs are active, and grow, at earlier epochs than lower mass black holes, an effect also seen 
bv lVestergaard fe Osmeil (|2009D and iSteinhardt fe Elvisl (|2010al) . 

4.3. The Broad Line Quasar Black Hole Mass Density, and Constraints on their Duty Cycle and Average Lifetime 

In Figure[5]we show the comoving mass density of SMBHs that reside in BLQSOs, as a function of redshift, pqso{z). 
The peak in the cosmic mass density of BLQSO SMBHs occurs at z ~ 2. Our constraint on the location of the peak 
in Pqso ( z) is consistent wi th the peak in th e mass density derived from the Large Bright Quasar Survey and high-z 
sample of lFan et al.l (|2001a| ). as calculated bv lVestergaard fe Osmerl (|2009l ). Previous work has not found any evidence 
for evolution in the Eddington ratio distribu t ion at z > 1 for the most massive BLQSOs (e.g., iMcLure fe Dunlopl 
[200l lVestergaardl[200l iKollmeier etaD [2001 1 Vestergaard fe Osmeri l2009). If there is no significant evolution in the 
Eddington ratio distribution at z > 1 for these systems, then we would expect the peak in their luminosity density 
to coincide with the peak in their bla ck hole mass density. The luminosity density of quasars peaks at z ~ 2 (e.g., 
Wol feT^[200l [Hopk ins et a l. 2007b), matching the peak in black hole mass density observed in our work. 

We can use the quasar BHMF to place constraints on the fraction of black holes that are seen in the BLQSO phase 
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Fig. 6. — Evolution in the mass density of SMBHs seen as BLQSOs, symbol s as in Figure [3] For comparison, recent estimates of the 
local mass density of all SMBHs suggest pbh(z = 0)~4x lO 5 M Mpc-3 fe.g. lYu fc LullSooa V The mass density of SMBHs in BLQSOs 
peaks at z ~ 2. 



as a function of Mbh', this quantity is commonly referred to as the BLQSO 'duty cycle'. Denote the duty cycle as 
6(Mbh,z). Then 

8{M BBt z) = * so % B "'*l (12) 
(/)bh{Mbh,z) 

Here, 4>bh(Mbh , z) is the BHMF of all SMBHs at a given redshift. Ignoring mergers of SMBHs, we can compute 
a lower-limit to the duty cycle by comparing the BHMF at a certain redshift with the local SMBH number density, 
since 4>bh(Mbh, z) < 4>bh(Mbh,0)- In Figure [7] we compute the lower limit of the BLQSO duty cycle at z = 1 
for SMBHs with M BH > 5 x 10 8 M Q . The duty cycle at z = 1 is constrained to be 6 > 0.01 at M B h ~ 1O 9 M , 
falling steeply to 5 > 10" 5 at Mbh ~ 10 10 M Q . The decrease in the duty cycle with increasing black hole mass may 
be another reflection of SMBH downsizing, with the most massive BLQSO SMBHs being active at earlier cosmic 
epochs. Alternatively, it may be due to the fact that the most massive SMBHs sp end a shorter amount of time in 
the broad line phase, as expected from some simulations of black hole feedback (e.g., Hopkins et al. 2006a). However, 
we note that the local number density of the most massive SMBHs is poorly constrained and subject to considerable 
systematic uncertainty (|Lauer et al.l 120071 ). and therefor duty cycle of BLQSOs with Mbh ~ 1O 1O M may be subject 
to considerable systematic error. Indeed, the observed fall-off in the lower limit on the duty cycle for the most massive 
systems may simply be due to the fact that we are overestimat ing the numbe r density of local SMBHs. 

The quasar lifetime is not well constrained empirically (e.g.. iMartinil 12004 ) . but we can use our estimated BLQSO 
BHMF to estimate the BLQSO lifetime. For simplicity, we assume that all BLQSOs of a given mass have a single 
lifetime, tsL- Furthermore, if a SMBH undergoes multiple episodes of BLQSO activity, then we assume that tsL is the 
same for each episode. These assumptions are unlikely to be true, but for the purposes of our work we may still think 
oUbl as a 'typical' BLQSO lifetime. Because the duty cycle is the probability of observing a SMBH as a BLQSO 
at a given Mbh and redshift we can relate tsh to S(Mbh, z). The probability of observing a SMBH as a BLQSO at 
redshift z is the product of the probability that a SMBH becomes a BLQSO along our line of sight at some point in its 
evolution with the probability that that SMBH is a BLQSO between cosmic ages t(z) — tsL and t(z), normalized by 
the probability that that SMBH is a BLQSO at t < t(z). The normalization results from the additional assumption 
that if a SMBH of mass Mbh exists at t(z), and if it goes through at least one BLQSO episode as some point in 
its growth, then it had to have gone through at least one BLQSO episode at t < t{z) in order to grow to Mbh by 
t(z). This is a reasonable assumption, at least for the most massive system, since all SMBHs at t{z) had to grow from 
much lower mass seeds, and if BLQSO activity occurs for a SMBH at some point in its life, BLQSO activity would 
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Fig. 7. — Lower bound on the duty cycle for BLQSO activity at z = 1 as a function of Mbh, symbols as in Figure [3] The duty cycle 
for Mbh ~ 10 9 M Q BLQSO SMBHs at z ~ 1 is 8 > 0.01, falling to S > 10" 5 for M BH ~ 10 10 M Q . The decrease in the duty cycle with 
increasing black hole mass is likely due to either cosmic downsizing or a shorter BLQSO phase for the most massive SMBHs. 



have occurred during this growth period. In other words, the assumption implicit in the normalization correction is 
that SMBH seeds of mass, say, M BB ~ 10 9 Af Q , do not simply emerge and then undergo BLQSO activity at a cosmic 
epoch later than t(z). The practical result of this assumption is that the BLQSO duty cycle of, say, SMBHs with 
Mbh ~ 10 9 Mq, can decrease from <5^1at,z^7to(5^ 10 -3 at z ~ 1. 

Denote the time that a BLQSO phase is initiated in a SMBH as r. Then, the probability that a SMBH BLQSO is 
seen between t{z) — tsL and t(z) is the same as the probability that the time that a BLQSO phase was initiated for 
that source occured at t(z) — t B i < t < t(z). Note that this allows for multiple BLQSO episodes, as multiple phases 
of BLQSO episodes simply alter the probability of observing a SMBH as a BLQSO at a certain redshift. The duty 
cycle is thus related to t B ^ as 



S(M BH ,z) = Pr(BLQSO\M BH ) 



Jt$_ tBL p(T\M BH , BLQSO) dr 



f* {z) p(T\M BH , BLQSO) dr 



(13) 



where, Pr(BLQSO\M B n) is the probability that a SMBH with a mass of Mbh is a BLQSO at some point in its 
evolution, and p{t\M B h, BLQSO) is the probability distribution of BLQSO episode initiation time for a given black 
hole mass. Note that Pt(BLQSO\Mbh) is not the probability that an object is currently seen as a BLQSO, but is 
the probability that it appears as a BLQSO to an observer on Earth at some point in its life. As mentioned above, 
p(t\M B h, BLQSO) is sufficiently general to include both multiple and single episodes of BLQSO activity, although 
the actual form of p(t\Mbh , BLQSO) will depend on the distribution of the number of BLQSO episodes a SMBH 
goes through. All values of Mbh in Equation (TT3")) refer to the mass of the SMBH at redshift z. 

If the distribution of r does not evolve significantly over a time t bl , then the integral in the numerator of Equation 
(|13p can be approximated as p{t\Mbh, BLQSO)tsL- Likewise, if p{t\Mbh, BLQSO) is approximately constant over 
tsL, then the distribution of BLQSO initiation times will be similar to the distribution of BLQSOs at the time we 
observed them. This is a reasonable assumption so long as tsL is short compared to the timescale for a significant 
change in the process that initiates BLQSO activity (e.g., mergers and other fueling events). Making this assumption, 
and rearranging Equation (I13|) . it follows that 



6{M BH ,z) J* iz) p(t'\M B H, BLQSO) df 
p{t{z)\M BH ,BLQSO)Pr{BLQSO\M B H) 



t B L 



(14) 
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Fi g. 8 . — Probability distribution, given the observed data and assumptions of our analysis, of the lower bound on the BLQSO lifetime 
fEq,[14|) calculated at z = 1 for Mbh = lO 9 Af0. The upper axis gives the lifetime in terms of the Salpeter timescale, which is tg a i p = 
4.3 X 10 7 yrs for a radiative efficiency of e r = 0.1. The lower-bound on the BLQSO lifetime is equal to the BLQSO lifetime if (1) all 
SMBHs of Mbh = 1O 9 M0 go through a BLQSO phase and (2) if the number densit y of these SMBHs at z = 1 is not significantly less 
than the local number density. There is evidence th at the latter condition is tr ue (e.g.. [Me rloni & Heinz 2008), while the former condition 
is expected if the SMBH's growth is self-regulated (Hopkins & Hcrnquist 2006). 

where p(t(z)\MBH, BLQSO) is the probability distribution of BLQSOs as a function of cosmic age, t(z). The term 
p(t\M BH , BLQSO) is related to the BHMF for BLQSOs according to 

<l>(M BH ,z) 1 
J™4>(M BH ,z) dM BH \ ' 

From Equation (|14p it follows that for the special case where all SMBHs experience a BLQSO phase, and where 
BLQSOs initiation times are uniformly distributed over the age of the Universe, then t B L = $(M BB , z)t(z). This 
is the definition of a quasar 'lifetime' commonly found in the literature; however, as the assumption of a uniform 
distribution of BLQSO initiation times i s inconsistent with the BL QSO BHMF or luminosity function, t B ^ will in 
general not equal S(M BB , z)t(z) (see also lHopkins fc Hern auist 2009[). 

Because we have assumed that t B L is the same for all BLQSOs of a given mass, and therefore must be the same at 
all redshifts, Equation (|14p may be calculated at any redshift. However, in reality we can only estimate a lower limit 
to t B L- This is because we do not know the fraction of SMBHs that will experience a BLQSO phase, although if the 
SMBHs growth is self-regulated we might expect Pr(BLQSO\M BB ) ~ 1. However, if some SMBHs of mass Mbh 
can never be seen as BLQSOs, possibly due to orientation-dependent obscuration, then Pr(BLQSO\M BB ) < 1. In 
addition, as discussed above, we can only calculate a lower limit to S(M BB ,z) by comparing the BLQSO BHMF at 
redshift z with the local BHMF of all SMBHs. Moreover, implicit in these calculations is the assumption that none of 
these quantities depend strongly enough on mass to vary significantly during the growth that occurs in the BLQSO 
phase. In Figure[5]we show the probability distribution of the BLQSO age of SMBHs with M BB = 1O 9 M , calculated 
from Equation (fl4|) at z = 1. We calculate t B L at M BB = 10 9 M Q because our sample is reasonably complete at this 
mass, especially at this redshift. In addition, we chose z = 1 because the local BHMF better approximates the z = 1 
BHMF than the BHMF at z > 1, therefore giving the tightest lower bound on the duty cycle at z = 1. We estimate 
t B i = 150 ± 15 Myr as a lower bound on the age of the BLQSO phase for SMBHs of Mp B = 10 9 Min. This estimate is 
roughly consistent with other est imates of quasar lifetimes (e.g. JYu fc Tremainell2002t iMartin i 2004; Goncalves et al.l 
2008| IHopkins fc Hernauisti[2009l) . 
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4.4. Constraints on The Most Massive Black Hole in the Universe 
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Maximum BLQSO M BH / M & 1 < z < 4.5 Redshift of Maximum BLQSO M BH 

Fig. 9. — Probability distribution of Mgu for the most massive SMBH that could be observed as a BLQSO at 1 < z < 4.5 (left), and the 
probability distribution for the redshift that this SMBH would be seen as a BLQSO (right). Here, we do not interpret M^^ x to represent a 
hard physical upper limit to Mgg, but rather it is the result of a finite number of black holes drawn from a mass function. The constraints 
on Mgjj® that we have obtained are consistent with previous observational work, but this is the first time that rigorous constraints on 
M^g x and its redshift have been obtained. The most likely value for Mgg x is M B H ~ 2.5 X lO lo Af , and this SMBH would likely be 
seen as a BLQSO at z > 2. 



The probability distribution for the maximum mass of a SMBH in a BLQSO may place important constraints on 
models of SMBH growth, and may be calculated directly from the BHMF. In the context of our work, we do not 
consider the maximum mass of a SMBH to be caused by a hard upper limit, above which it is impossible to make 
a more massive black hole, but rather the result of a finite number of black holes drawn from a mass function. The 
probability that the maximum mass of a SMBH in a sample of N BLQSOs is less than M. is simply given by the 
probability that all N SMBHs have M BB < M: 

Pr(M^ x <M\N) = [Pr(M BH < M)f (16) 

Note that N is the total number of BLQSOs that could be observed in an all-sky survey with no flux limit; i.e., TV is 
the normalization of the BLQSO BHMF. The term Pr{M BH < M) is calculated from the BHMF as 

i r M r°° fdV\ 
Pr{M BH < M) = — J J <f>{M BH ,z)(—J dzdM BH . (17) 

The probability distribution of M B jj x is then found by differentiating Equation (TTB1 with respect to M. and evaluating 
the result at M = M^ x : 

p(M^ x \N) = [Pr(M BH < M BB X )] J™ <t>(M BH , z) [^\ dz. (18) 

The posterior probability distribution for Mg^ x , given the observed data, can be calculated by averaging Equation 
(fT8]) over the MCMC output. 

In Figure |9] we show the posterior probability distribution of the maximum SMBH in a BLQSO at 1 < z < 4.5. 
Here, Mjj^ should be interpreted as the maximum mass of a SMBH in a BLQSO that would be observed in an all-sky 
survey without a flux limit over the redshift interval 1 < z < 4.5. Thus, Mg^ x is a lower bound on the mass of the 
most massive SMBH in the Universe. In addition, in Figure [9] we also show the probability distribution of the redshift 
at which the quasar with Mjffi* would be found. As can be seen from these figures, the maximum mass of a SMBH 
in a BLQSO at 1 < z < 4.5 is M B ^ X - 3 x 10 10 M Q . The probability distribution for the redshift of M^ x is rather 
broad, but we constrain the redshift for the BLQSO hosting this SMBH to be z > 2. 

Our results are in agreement with what others have fou nd using samples of broad li ne mass estimates (e.g., 
lVestergaard1l2004t iVestergaard et al.H2008b iNetzer et alJl2007l) althou gh lLabita et all (|2009al lbl) find a smaller value of 
Mf™ ~5x 10 9 M©. However, in the model of lLabita et a l. (2009a) Mj,^p is a parameter determining the shape of 
the mass function, and is not the actual maximum black hole mass o f a sample of objects drawn from the distribution 
of black hole mass; i.e., there is nothing in the statistical model of lLabita et al . (2009a) that prevents objects with 
M B h > Mg£ x . As such, the actual realized maximum mass in a large sample of objects will be greater than the 
value of estimate d using the model of l Labita et al.l (|2009al ). as we have found. Considering this, our results are 

consistent with those of lLabita et al . (2009a.^ 

These highly massive black holes represent the extremes of black hole growth, and thus are important for constraining 
models of black hole growth. Furthermore, these massive black holes offer the best chance of probing the faint end of 
the BLQSO Eddington ratio distribution. Therefore, it is of use to investigate how large and deep a survey must be 
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Fig. 10. — Expected number counts of BLQSOs at 1 < z < 4.5 with black hole masses 5 X 10 s < Mbh/Mq < 5 X 10 9 (solid line) and 
5 X 10 9 < M BH /M Q < 5 X 10 10 (dashed line) in a 1 deg survey as a function of limiting magnitude. One would need to search an angular 
area Q > 10 deg 2 in order to expect to find at least one BLQSO with Mgjj ~ 10 10 Mq. In addition, one does not become complete at 
^Ibh ~ 5 x 10 8 M Q until i ~ 23. 

in order to find an useful number of these objects. In Figure [TU] we show the expected number of high mass SMBHs 
in BLQSOs at 1 < z < 4.5 in a 1 deg 2 survey as a function of limiting i-magnitude, assuming our best-fit statistical 
model. We stress that these are the expected number counts for the true black hole masses, and not the mass estimates. 
Because the errors in the mass estimates will scatter more sources into higher mass bins than into lower mass bins, 
the number of sources with estimated mass in a mass bin will be larger than the actual number of sources. 

Any survey with a flux limit of i < 19 should be able to detect BLQSOs with Mbh ^ 5 x 1O 9 M ; however, the 
survey must have an area of Q, > 10 deg 2 to expect to detect at least one of these objects. Similarly, in order to 
get a large number of BLQSOs with Mbh }t 5 x 10 8 M Q , one needs a deep, large area spectroscopic survey. For 
example, the BLQS O spectroscopic samples from th e COSMOS survey ([Scoville et aLll2007f ) cover « 2 deg 2 at i < 24 
(|Trump et al.l 120071) and i < 22.5 for zCOSMOS (jMerloni et al.ll2010D . and, ignoring redshift incompleteness, are 
therefore expected to contain ~ 60 BLQSOs at 1 < z < 4.5 with masses 5 x 10 s < Mbh/Mq < 5 x 10 9 , but none with 
masses Mbh > 5 x 1O 9 M . 



As discussed in § 13. 1[ we also estimate the distribution of the ratio of the Eddington ratio to the bolometric correction 
at 1350A, TEdd/Cisso- Under the assumption that T Edd/ C\^o follows a mixtu re of log- normal distributions , with mean 
values that depends linearly on log Mbh, and that C1350 = 4.3 for all sources (|Vestergaard fc O smer 2009), our model 
implies that the geometric mean Eddington ratio increases with Mbh according to: 



The dispersion in L/LEdd decreases with increasing Mbh, having a value of ~ 0.4 dex at Mbh ~ 10 8 M Q , and 
decreasing to ~ 0.3 dex at Mbh ~ 10 9 M Q . Both the dispersion and slope are smaller than what was found by KVF09, 
who analyzed a sample of BLQSOs from the Bright Quasar Survey at z < 0.5, suggesting possible evolution in the 
Eddington ratio distribution. However, the statistical significance of a difference in the slopes is marginal at best, as 
the differences are only significant at ~ 2a. 
In Figure[TT]we show the implied distribution of BLQSO Eddington ratios assuming a constant bolometric correction 



4.5. Implied Eddington Ratio Distributions 




(19) 
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Fig. 11. — Inferred Eddington ratio distribution, assuming the mixture of log-normals form of Equation J2I and a bolometric correction 
of C1350 = 4.3. The vertical line marks the 10% completeness limit for the SDSS DR3 sample at z = 1; Figure l4l shows that the 
10% completeness limits are similar for z ~ 1,3, and 4, but shallower for z ~ 2. Our inferred Eddington ratio distribution peaks near 
L/LEdd ~ 0.05 and has dispersion of ~ 0.4 dex, although the peak fall below the 10% completeness limit and should be interpreted with 
caution. Our inferred Eddington ratio distribution is shifted toward lower values of L/L^dd an d nas a higher dispersion than previous 
work, which did not correct for incompleteness. Our estimated Eddington ratio distribution suggests that most SMBHs in BLQSOs are 
not radiating near their Eddington limit. 



of C1350 = 4.3 ([Vestergaard fc Osmeii [20091. As can be seen, the estimated distribution of quasar Eddington ratios 
peaks at L/LEdd ~ 0.05. However we note that the location of the Eddington ratio peak falls below the 10% 
completeness limit of the SDSS, and thus the exact location of the peak is highly uncertain. Our inferred Eddington 
ratio distribution is broader and shifted towards smaller values of L/L Edd than what has been foun d in previous work, 
which were not able to fully correct for incomp leteness in black hole mass and Edding ton ratio (e.g..lMcLure fc Dunlopl 
12004 lVestergaardll200llKollmeier et al.ll2006h . However, it is consistent wit h what IShen et al I 1)20081) found using a 
similar approach to ours. The major differences between our work and that of iShen et al.l ( 20081) is that they assume a 
power-law distribution of Mbh, and they constrain the BHMF and Eddington ratio distribu tion by visu a lly ma tching 
the model distributions to the observed distributions. Our results, as well as the results of IShen et al . (2008), show 
that the narrow distribution in quasar Eddington ratios seen in previous work, and peaking at L/LEdd ~ 0.25, is due 
to uncorrected incompleteness. Indeed , deeper samples of BLQSOs have observe d Eddington ratio distributions that 
are broader and peak at lower values ()Gavignaud et al.l 120081 : [Trump et al.ll2009D . Therefore, we conclude that most 
BLQSOs are not radiating at or near the Eddington limit, and that there is a large dispersion in Eddington ratio for 
BLQSOs. 

Rec ently, ISteinhardt fe Elvisl (|2010aj ) have analyzed the SDSS DR5 sample of mass estimates calculated by lShen et al.1 
(2008), and concluded that there is a dirth of objects at high Eddington ratio, and that there is a redshift-dependent 
systematic decrease in Eddington ratio with increasing black hole mass for BLQSOs in the high L/LEdd tail of the 
distribution. They have called this effect a sub-Eddington boundary. While we also find evidence that BLQSOs 
radiating near the Eddington limit are rare, as did KVF09, the dependence on black hole mass is seemingly in contrast 
to what we find here, in that we find that the average Eddington ratio increases with Mbh, a ssuming a constant 
bolome teric correction. The most important difference between their work and ours is the fact that ISteinhardt Sz Elvisl 
(2010 a|) analyze seperate re dshift bins, while we do not model a redshift dependence in the Eddington ratio distribution. 
Stcinhard t fc Elvisl (|2010a|) did not perform any calculations for the whole sample, averaged over all redshifts, making 
a more direct comparision difficult, and thus it is unclear how discrepant our conclusions are. 

Another potential contributor to this discrepancy is the handling of the statistical error in the mass estimates. 
Statistical errors (e.g., measurement error ) have the effect of flattening slopes and correlations when one analyzes the 
quantity contaminated by the error (e.g.. lAkritas &; Bershadvl 119961 ; lKellvl l2007; Kel ly fc Bech told 20Q3), due to the 
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fact that the distribution of the estimated quantity is a biased estimate of the distribution of the quantity of interest. 
Because the mass estimates are contaminated by statistical error, the dependence of luminosity on the mass estimates 
will be natter than the intrinsic dependence of luminosity on Mbh- In particular, the statistical error will cause some 
masses to be underestimated, and some to be overestimated. The objects which appear to have the highest masses 
at a given luminosity will have the highest overestimates, and thus will have their Eddington ratios underestimated. 
As a result, there will be an apparent dirth of high Eddington ratio objects at the highest masses. Similarly, there 
will be an apparent excess of high Eddington ratio objects in the low mass bins. Our method estimates the intrinsic 
dependence of luminosity on Mbh by, in a sense, 'deconvolving' the distribution of the mass estimates and luminosity 
with the distribution of the error in the mass estimates; this is also why we find a somewhat steeper decline in the 
mass function at the highest masses, as compared to the estimate reported by [Vcstcrgaard et al. (20Q3). 

Although the statistical errors in the mass estimates can have the affect described above, possibly contributing to 
the different conclusions, it is unlikely that the statistical errors alone are sufficient to account for the discrepancy, 
and other possibilities include differences in the methods and scaling relationship employed for estimating the masses, 
and other differences in the methods of data analysis employed. However, further investigation is beyond the scope of 
our work. In add ition, we note that our sam ples only overlap at 1 < z < 2^1 . and we cannot compare with the low 
redshift results of Stcinhar dt fc Elvisl (|2010 a). where they find the greatest evidence for these effects. 

Although we have used a flexible form for p{L\Mbh), our estimated distribution for L/LEdd may be affected by 
significant systematics, as it relies on the assumption of a constant bolometric correction. There is ey idence that the 
bolom etric correction de pends on both Eddington ratio (jVasudevan fc Fa bian 20Q2 lYoung et al.ll2010l ) and black hole 
mass (|Kellv et al.l 12008). and the increase of L /LEdd with Mbh we find may instead be a reflection of a decrease in 
the bolometeric correction with Mbh- Indeed, iKellv et al] (|2008l ) find that the ratio of optical/UV flux to X-ray flux 
increases with Mbh, implying that the bolometeric correction to the UV flux decreases with Mbh, and therefore we 
would infer an increase in L/LEdd with Mbh assuming a constant bolometeric correction, as we do. In addition to 
these issues, we also note that the SDSS is at least partially incomplete at all Eddington ratios at z > 1 El as shown 
in Figure SJ and further work is needed using deeper surveys to confirm these results. 

5. DISCUSSION 
5.1. Connection with Quasar Lightcurves 

A significant amount of recent work has suggested that quasar feedback regulates the growth of SMBHs and affects 
the large-scale evolution of its host galaxy. Understanding the quasar lightcurve is thus of fundamental importance 
for understanding black hole growth and feedback, as well as placing constraints on the physics of quasar accretion 
flows. The distribution of luminosity at a given black hole mass can be related to the quasar lightcurve, and therefore 
one can use the estimated distribution of luminosity at a given black hole mass to place some empirical constraints 
on models for quasar lightcurves. We address this in the following section, and argue that our estimated Eddington 
ratio distribution is consistent with models where the BLQSO phase represents the final stages of a SMBH's growth, 
in which the QSO lightcurve is expected to decay until the object no longer appears as a BLQSO. 

For a given fueling episode, denote the quasar lightcurve as L(t), and the mass fueling rate to the SMBH as a 
function of time as M(t). The fueling rate, M(t), gives the rate at which matter is externally supplied to the BLQSO 
accretion disk. The probability distribution of quasar luminosity at a given Mbh, M, and t is p(L\Mbh, M, t), and the 
probability distribution of the fueling rate at a given Mbh and t is p(M\Mbh ,t). The distribution p(L\Mbh, M,t) 
is determined by the physics of the accretion disk, since the quasar luminosity is produced by viscous stresses in 
the accretion disk. For example, iSiemiginowska &; Elvisl (|1997l ) calculate p{L\MsH,M,t) implied by the model of 
iSiemiginowska et~a l. (1996) for a thermal- viscous accretion disk stability. The distribution p(M\Mbh, t), on the other 
hand, depends on the physics and stochastic nature of the fueling mechanism and quasar feedback. 

Because the quasar luminosity is determined by the physics of the accretion disk, it is reasonable to assume that 
given Mbh and M , L is independent of time. This does not imply that the quasar lightcurve does not vary with time, 
as it certainly does, but rather that knowing the value of t does not give us any additional information on the value 
of L when we already know Mbh and M at a given t. We therefore drop the explicit dependence of p(L\Mbh, M,t) 
on time. 

The distribution of luminosity at a given Mbh is calculated as 

rt +tBL 



P (L\Mbh) = / p{L\M,M BH ) 



p{M\t,M B H)p{t\M B H) dt 



to 



dM, (20) 



where p(t\MBn) is the probability of seeing a BLQSO at time t, given Mbh, and tsL is the length of time that a quasar 
would be classified as a BLQSO. Here we have defined the broad line phase of quasar activity to start at t — to. The 
term p(t\MsH) is related to the model for black hole growth, since p(t\MsH) oc p(Mbh \t)p(t). Because we observe 
quasars randomly during their BLQSO phase, p(t) is uniform and p(t\MsH) °c p(Mbh \t)- In general, BLQSOs with 
larger values of Mbh are more likely to be seen later in their lifetime, i.e., at larger values of t. 

1(1 Stcinhardt k Elvis (2010a) also included quasars at z > 2 in some of their plots, but their calculations and conclusions were based on 
quasars at z < 2 because of concerns regarding the validity of using C IV in the mass estimates. 

11 As discusses earlier, this is merely a reflection of the fact that the mass function increases steeply toward lower masses, which are hard 
to detect even if these low mass SMBHs are radiating near the Eddington limit 
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As an alternative to Equation (f2T))). p(L\Mbh) may be directly calculated from the amount of time that a BLQSO 
spends above a given luminosity, as a function of Mbh- I n this case, the term v( L\Mbh) is directly related to the 
'luminosity-dependent' lifetime interpretation advocated bv lHopkins et a l. (2005b. c), where the quasar 'lifetime' is the 
amount of time that a quasar spends above a given luminosity. Assume that p(t\MsH) oc 1 for BLQSOs, and denote 
T{L\Mbh) to be the amount of time a BLQSO spends above a luminosity L, as a function of Mbh- Then 

dT(L\M BH ) 



p(L\M BH ) oc 



dL 



(21) 



For example. iHopkins fc Hernquistl (|2009D suggest that p(L\Mbh) can be well approximated as a Schechter function, 
which generalizes the distributions implied by several common models of quasar lightcurves often found in the literature. 

5.1.1. Effects of Evolution in the Quasar Fueling Rate 

Recent work on quasar fueling and feedback suggests that the BLQSO phase occurs at the end of the fueling event, 
during which feedback energy from the AGN is able to unbind the surrounding gas, removing the obscuring material 
and halting the black hole's growth. Within the context of this model, the SMBH does not experience a large fractional 
increase in mass beyond the value of Mbh that is observed during its time as a BLQSO. Therefore, we approximate 
p(t\MBn) as being uniform and independent of Mbh, i-e., p{t\MsH) = {/tBL- In addition, this model predicts that 
during this so-called 'blow-out' phase, the fuel supply is set by the evolution of the blastwave caused by the injection 
of ene rgy from the SMBH, and is expected to be of the form M(t) oc i _/3 ([Hopkins & Hcrnquist 2006; Hopkins et al.l 
l2006a| ). Alternatively, if the fuel supply is suddenly removed, then the accretion rate during the broad line phase is 
due to evolution of the viscous accretion disk. The fueling rat e also has a p ower- law form under viscous evolution of 
the disk, but with a different value for /3 (|Cannizzo et al.lll990HPringlel fl991). Under these models, it is reasonable to 
approximate the fuel supply as a deterministic process, M(t) oc . If the evolution of the fuel supply is a deterministic 
and one-to-one process, and assuming that p{t\MBn) °c I/^blj then p{M\Mbh) is 

<M\M \ 1 dt ( M ' M BH) , * 

P{M\Mbh) = -T- , (22) 

tBL dM 

where t{M ,Mbh) is the time implied by M(t, Mbh), i-e., t(M, Mbh) is the inverse function of M(t, Mbh)- If the 
quasar lightcurve is a many-to-qne fun ction, then the right side of Equation (|22l) is replaced by a sum of derivatives, 
as in Equation (15) of lYu fc Lul (|2004j ). Inserting Equation ([22]) into Equation (|20|) and dropping the time integral 
give£3 

p(L\M BH ) = t~ !°° p(L\M,M B H) dtiM :*! BH) dil - ( 23 ) 
tBL Jo dM 

For models where M(t) oc , a power-law distribution of fueling rates follows from Equation (|22|) . p(M\Mbh) oc 
Af-U+i//3) ( see a i so Eq. 43 in lYu &; Lull2008[ ). Therefore, assuming M(t) oc t~ p gives 

p(L\M B h)(x p(L\M,M B h)M- {1+1/0 ' ) dM. (24) 



JO 

The feedback-driven model of [Hopkins fc Hernqui st ( 2006) predicts a value of [3 ~ 2, while the disk evolution model 
predicts f3 ~ 1.2 (jCannizzo et alJll990t lYu et alJl2005 iKing fc Pringlell2007t ). implying that p(M\M BH ) oc A/" 1 - 5 or 
p(M\M B h) oc M- 1 - 67 , respectively. 

5.1.2. Effects of Time- Dependent Accretion Disks 

It is apparent from Equation ([24]) that the effect of time-dependent accretion disk behavior is to create a distribution 
of luminosities at a given Mbh that is flatter than that expected from a simple power-law decay in the fueling rate. 
While the BLQSO fueling rate may be approximated as a deterministic process, qu asar lightcurves are observed to 
exhibit aperiodic and stochastic variations across all wavelengths (for a review see lUlrich et al.lfl997h . Therefore, 
the quasar lightcurve will be stochastic, i.e., the value of the luminosity is not completely determined by Mbh, and 
M. From Equation (|24|) it can be observed that if we assume a deterministic relationship between L and Af, i.e., a 
time-steady accretion disk, as much previous work on SMBH growth and feedback has assumed, then L oc M and 
the quasar lightcurve simply traces the fueling rate evolution. In this case, p(L\Mbh,M) is a delta function and 
p(L\Mbh) oc _L~( 1 + 1 /' 3 ). However, the time-dependent and stochastic nature of the accretion disk emission broadens 
the observed luminosity function beyond that implied solely from the BLQSO fueling evolution. 

Although quasars are observed to vary at all wavelengths, we focus our rem aining discu s sion on optical variability, 
since we are interested in the distribution of optical luminosities in this work. IKellv et al.l (|2009aD found that quasar 
optical lightcurves on timescales < 7 yrs in the rest frame of the quasar are well-described by a Gaussian process 

12 This follows from p(L\M BH ) = f p(L\M , M BH )p(M\M BH ) dM. Alternatively, we could have arrived at Eq.[23] by taking 
p(M\Mgn,t) to be a delta function centered at the value of M determined by t, and directly doing the integration in Eg. 1201 . 
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on the logarithmic scale (see also iKozlowski et al.l l2010h . with characteristic timescale consistent with accretion disk 
thermal time scales. They suggested that quasar variability on these time scales is due to a turbulent magnetic field 
in the accretion disk, which drives changes in the radiation energy density of the disk, as also seen i n 3-dimensional 
magneto-hydrodynamic simulations of radiation-pressure dominated accretion disks (|Hirose et a l. 2008). Extrapolation 
of their best-fit lightcurves imply that flux variations for this particular process have a standard deviation ~ 0.1 dex 
on timescales long compared to the disk thermal time scale, independent of Mbh- If this is the only source of 
quasar optical variability, then p(L\Mbh, M) is a gaussian distribution centered at L opt = C opt e r Mc 2 with standard 
deviation ~ 0.1 dex. Here C opt is the bolometric correction to the optical luminosity (which likely depends on 

Mbh and M), e r is the radiative efficiency, and c is the speed of light. This process only produces a small amount of 
broadening in the luminosity distribution, relative to that implied solely from the BLQSO fueling evolution. Therefore, 
p(L\Mbh) oc L~( 1+1 /^ is a reasonable approximation when the fuel rate decl ines as a power-law , so long as there are 
no other variability components in the accretion disk beyond that observed bv lKelly et al.1 (|2009aT ) . However, this does 
not avoid the issue of a non-constant bolometric correction. 

Variability on longer time scales may be driven by acc retion disk instabilities. Accretion disks based on the standard 
a-prescription for viscosity (jShakura fc Sunvaevl I19T3T ) are subject to a number of accretion disk instabilities that 
can potentially have a significant effect on the quasar lightcurve. At high accretion rates (m > 0.025) a radiation- 
pressure instability may operate on time scales > 10 4 yrs dCzernv et al.l [20091). while at lower accretion rates a 
thermal-viscous io nization instability (e.g., iLin fc Shields! 119861 : iSiemiginowska et al.|[l996t iMenou fc Quataertl 120011 : 
Uaniuk et al.ll200"4l ) may operate on time scales > 10 6 yrs. The ionization i nstability has been invoked to explain the 
outbursts seen in dwarf novae and soft X-ray transients (for a review see !Lasotall200lD . and the r adiation pressure 
instability has been invo ked to expl ain the outbursts seen in the microquasar GRS 1915+105 fe.g- INavakshin et al.1 
2000; Uaniuk et al.ll2000T ). Moreover. iGoodman fc Tai3 (|2004D suggest that stars may also form in AGN accretion disks 
due to fragmentation in the outer edges of the disk. While it is apparent that accretion disks around galactic sources 
exhibit instabilities, it is currently unclear if and how these instabilities operate in AGN, as the time scales involved are 
too long to observe transitions between quiescence and outburst. Furthermore, theoretical predictions of lightcurves 
based on these instabili ties vary, and the existence and importance of disk stabil ities can depend on how the viscosity 
is parameterized (e.g., iStella fc Rosnerl 119841; iSzuszkiewiczl 119901 : iMerlonil [2003). how much of the accretion energy 
is dissipated in a hot corona (e.g.. iSvensson fc Zdziarskilll994f ). or if there is an outflow or an advection dominated 
accretion flow (ADAF) (e.g.. lHameurv et al.ll2009| ). 

The location in the disk where the instability occurs is important, as p(L\M, Mbh) is with respect to the optical flux 
in our model. Various instabilities are expected to operate in different regions of the disc, sometimes in a stochastic 
manner, and thus may or may not significantly affect the optical flux. The term p(L\M , Mbh) is also conditional 
on the source being seen as a BLQSO, and if the disk instabilities cause the ionizing continuum to disappear during 
quiescence, thus causing the broad emission lines to disappear as well, then p(L\M , Mbh) is only with regard to 
the distribution of optical luminosities during a disk outburst. Considering these issues, and the current uncertainty 
regarding the importance of disk instabilities, we consider it beyond scope of this paper to fully investigate the effects 
of time-dependent behavior in the accretion disk. However, in light of this di scussion, disk instabilities can potentially 
have a significant effect on the distribution of luminosity at a given Mbh (Sicmiginows ka fc Elvisl[l997f ). Thus, the 
distribution of luminosity at a given Mbh is likely altered beyond a simple power-law expected if the optical flux is 
trivially related to the external fueling rate. 

In spite of these issues related to the detailed physics of the accretion disk, our estimated p(L\Mbh) is qualitatively 
consistent with models which predict the BLQSO phase to occur while the quasar accretion rate is decaying, as we find 
that high Eddington ratio objects are rare, and that the number density of BLQSOs increases steeply toward lower 
values of L/LEdd- In addition, although our estimated Eddington ratio distribution continues to rise toward lower 
L/LEdd, our sample is too incomplete to rigorously probe the low Eddington ratio region of the distribution, and our 
estimated L / L,Edd distribution in this region is heavily dependent on our parameteric form used to extrapolate beyond 
L/Lsdd S 0-1- Thus, we cannot test if the Eddington ratio distribution continues to rise until L/LE dd ~ 10~ 2 , after 
which a steep decay in p( L / L Edd) might occur due to the disappearance of the broad emission lines fe.g. lChurazov et al.l 
2005: lTrump et al".ll2009T ). 



5.2. Implications for the Growth of Supermassive Black Holes 

Our results in this work are broadly in agreement with recent observational and theoretical work suggesting that 
SMBH growth and spheroid formation have a common origin. In particular, models where mer gers dominate black 
hole growth predict that major mergers of gas-rich galaxies initiate a burst of star formation (Mihos fc Hernquistl 
1994a| Il996f) . during which the black hole undergoes Eddington-limited obscured growth. Eventually the black hole 
becomes massive enough for radiation-dr iven feedback to unbin d the surrounding gas, halting the accretion flow and 
revealing the object as a BLQSO (e.g., iHopkins et al.l [200 6b). As the activity further declines, the remnant will 
redden and become quiescent, satisfying the b lack hole-host galaxy correlation and leaving a dense stellar remnant 
from the starburst (M ihos fc Hernqu ist 1994 b]). This dense stellar remnant is identifiable as a second compo nent in 
the lig ht profiles of elliptical galaxies ()Hopkins et al.l [2008b. 2009c. d|). Techniques based on the argument of iSoltanl 
(11982T) have concluded that the SMBH accr etion rate density of the Universe peaks at z ~ 2 fe.g.. lMarconi et al.ll2004t 
iMerloni "fc Heinz 2008; Sha nkar et al.|[2009T) . in agreement with the observed peak in the SMBH luminosity density of 
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the Universe at z ~ 2 fe.g. JWolf et al.ll2003l; IHopkins et al.ll2007bD . We have found that the peak in the BLQSO cosmic 
mass density also occurs at z ~ 2, in broad agreement with models where the SMBH undergoes Eddington-limited 
growth up to the BLQSO phase, at least on a cosmological level. 

In addition, our results suggest that a lower bound on the typical length of time for which a massive (Mbh ~ 10 9 M Q ) 
SMBH could be seen as a BLQSO during a single BLQSO episode is tg^ > 120 Myr, i.e., a few Salpeter times. If 
all M B h ~ 1O 9 M SMBHs go through a BLQSO phase, and if the number density of SMBHs with M BH ~ lO 9 Af 
does not significantly increase from z = 1 to z = 0, then Equation (|14l) is no longer a lower bound and tsh ~ 150 
Myr represents an estimate of the lifetime of a single BLQSO phase. Indeed, if the SMBH's growth is self-regulated 
then there is good reason t o assume that most SMBHs go through a BLQSO phase at the end of their growth 
(jHopkins fc Hernquistl [2006h . as the obscured/unobscured dicho tomy represen ts an evolutionary sequence, at least 
at z > 1. Furthermore, recent work employing the argument of iSoltanl (|1982h suggests that the number density of 
M B h ~ 10 9 Af© SMBHs does not significantly increase from z = 1 to z = (jMerloni fc Heind [20081 ) . These two 
considerations imply that tsL ~ 150 Myr should represent a reasonable estimate of the length of the BLQSO phase. 

Our estimat e d BLQ SO lifetime of tsL ~ 150 Myr is longer than the value of tsL ~ 10-20 Myr predicted by 
Hopk ins et al.l (|2006bl ). Part of this discrepancy may be due to the fact that they defined a BLQSO as having a 
B-band luminosity brighter than some fraction of the host galaxy's, while we define a BLQSO as simply having broad 
e mission lines along our line of sight. Therefore, if a BLQSO has a decaying lightcurve, as assumed in the model 
of IHopkins et al.l (|2006bD . then by our definition a SMBH may still be considered a BLQSO even after the B-band 
luminosity has fallen below its host's. This would imply that the predicted BLQSO lifetime under the definition of 
IHopkins et al.1 (|2006bl ) would be shorter than our estimated value, as is the case. However, i t is unclear if the lifetime of 
the BLQSO should be a factor of ~ 10 shorter under the definition of lHopkins et al.l (|2006b| ). Unfortunately, we cannot 
invert our estimated Eddington ratio distribution to a quasar lightcurve to test this, as the inversion is not unique, 
and we would have to assume a distribution o f host g alaxy B-band luminosity during the BLQSO phase at z ~ 1, 
which is poorly constrained; as Hopkins et al.l (|2006bf) point out, uncertainty in the ratio of host galaxy to quasar 
B-band luminosity also introduces considerable systematic uncertainty in their predicted BLQSO lifetime. Moreover, 
our estimated lifetime is subject to the assumptions outlined in § 14. 3[ which may introduce additional systematic 
uncertainties of a factor of a few. Co nsidering this, a more quantitative comparison between our estimate and the 
prediction from Hop kins et al.l (|2006bf) is difficult. 

Our estimated BLQSO BHMF suggests a typical value of Mbh ~ 10 8 Mq for BLQSOs. Theoretical estimates of 
the mass distribution of SMBH se eds suggets that typica l values for the seed mass should be M BH ~ 10 5 -10 6 M Q at 
z ~ 10 (lLodato fc Nataraianl[2007t iPelupessv et~aT1l2007t iVolonteri et~aTll2008l ). and grow to M BH ~ 10 8 M Q by z ~ 3 
(|Volonteri fc Nataraianl 120091 ). consistent with our results. If the black hole's growth is Eddington-limited, then at 
least several Salpeter times are required to grow a seed black hole to M B h ~ 10 9 M Q from a M BH ~ 10 6 M Q seed, 
which is inconsistent with our estimated BLQSO lifetime. Furthermore, we find that most SMBHs in BLQSOs are not 
radiating near the Eddington limit, with a typical value being L/LEdd ~ 0.1 at Mbh ~ lO 9 Af0, as we might expect 
if BLQSOs are seen during a phase with a decaying fueling rate. Because most BLQSOs are radiating at considerably 
less than the Eddington limit, this therefore suggests that a factor of ~ 10 longer is needed to grow these SMBHs from 
seeds of Mbh ~ 1O 6 M0, i.e., a growth time scale of ~ 70 Salpeter times. This corresponds to the age of the universe 
at z ~ 2, implying that sources observed at z > 2 had to have been accreting at higher rates earlier in their growth. 

An inferred growth time scale of ~ 70 Salpeter times is significantly longer than our estimated BLQSO lifetime 
at z ~ 1, implying one of a few explanations. First, it implies that if we assume that BLQSO SMBHs spend all of 
their growth as a BLQSO, then our assumption that all SMBHs go through a BLQSO phase along our line of sight 
is incorrect. In other words, this implies that some of the SMBH population that is growing at any redshift could 
never be observed by us to be a BLQSO at any time. Because the lifetime is related to what we can calculate from 
the BHMF by Equation ([13]). and if the SMBH spends all of its growth as a BLQSO, then an inferred SMBH growth 
time scale of t B L ~ 70t sa i p implies that only ~ 4% of SMBHs ever go through a phase where we could observe them 
as BLQSOs at any time. This is an unrealistically low number, and thus we conclude that tsL is shorter than the 
growth time. Furthermore, it also illustrates that even if we can never observe a large fraction of active SMBHs, say 
~ 50 — 75%, possibly due to orientation-dependent obscuration, then our estimated tsL is only underestimated by a 
factor of ~ 2-4. 

Another possibility is that SMBHs in BLQSOs at z ~ 1 built up their mass via numerous fueling events of length 
tsL ft 150 Myr. In this case one could grow a SMBH from a seed mass of Mbh ~ 10 6 Af Q at z ~ 10, say, to a mass 
of Mbh ~ 1O 9 M by z ~ 1, without the need for an earlier phase of obscured and accelerated growth. Similarly, we 
may have incorrectly assumed that the lifetime of the BLQSO phase is short compared to the time scale needed for 
any significant change in the BLQSO triggering time distribution, as this assumption was needed in order to use the 
observed distribution of BLQSOs as an estimate of their triggering distribution. This assumption may be incorrect if 
BLQSOs undergo a single long growth phase from z ~ 10 to z ~ 1. However, this possibility only exists for z < 1.7, 
since at higher redshifts the growth time for objects radiating at L/LEdd ~ 0.1 is longer than the lookback time to 
z ~ 10, and thus BLQSOs that are observed at z > 1.7 would have had to experience an earlier phase of obscured 
growth at an enhanced accretion rate. Moreover, if SMBHs that are seen as BLQSOs at z ~ 1 with Mbh ~ lO 9 Af 
are grown in a single long BLQSO episode, or numerous repeated ones of tsL ~ 150 Myr, this would also imply 
that BLQSOs that are seen at z < 1.7 would have had a different fueling mechanism than BLQSOs that are seen at 
z > 1.7, and it is unclear why this should be true. Therefore, we conclude that while part of the mass of SMBHs in 



22 



Kelly et al. 



BLQSOs may have been accumulated via multiple BLQSO episodes, it is unlikely that most of these SMBHs did not 
also experience a phase of obscured growth at an enhanced accretion rate. 

If most SMBHs go through a BLQSO phase along our line of sight at some point in their growth, the fact that our 
estimated BLQSO lifetime is short compared to the SMBH growth time implies that SMBHs spend a significant amount 
of their time growing in a non-BLQSO phase, such as an obscured phase. Note that this argument is unaffected by the 
fact that at any z we miss a considerable fraction of the SMBH population that is growing (e.g., due to obscuration), 
so long as these missed SMBHs undergo a BLQSO episode at some point in the ir life. The conclusion that a significant 
amount of growth occurs in an earlier obscured phase was also reached by iTreister et al.l (|2010D . and is expected 
from self-regulation models models where mergers or other triggering mechanisms fuel and initiate Eddington-limited 
accretion and obscured SMBH growth, until the SMBH becomes massive enough to unbind the ambient gas, revealing 
it as a BLQSO for t B L ~ 150 Myr. Within this interpretation, the BLQSO black hole mass density shown in Figure E] 
is proportional to the rate at which the relic SMBH mass increases as a function of redshift. Furthermore, considering 
that the growth time scale for a SMBH radiating at L/LEdd ~ 0.1 to grow from 10 6 M Q to 10 9 M Q corresponds to the 
age of the universe at z ~ 2, we also conclude that SMBHs accrete at a significantly higher rate during the earlier 
obscured phase, as compared to the BLQSO phase. 

In this work we have found that the most massive SMBH that could be seen as a BLQSO is M^ x ~ 3 x 1O 1O M0 
and would most likely be observed at z > 2. These constraints on Mg£ x are consistent with recent cosmological 
simulations of SMBH growth, as well as expectations from black hole feedback models. Cosmological simulations that 
follow the growth of SMBHs in brig ht z ~ 6 quasars have been able t o grow SMBHs to Mb h ~ 10 9 M Q by z ~ 4 
(|Li et alJl2007tlDi Matteo et al.ll2008h and M BH - 2 x 1O 1O M by z ~ 2 (jSijacki et al.ll2009f ). Similarly, considerations 
based on quasar feedback and self-regulated SMBH growth, combined with the local distribution of bulge velocity 
dispersion, also suggest a value of M§fi x - 10 10 M Q (jNataraian fc Treisterll2009f ). 

If the SMBHs growth is self-regulated, then the final mass of the SMBH a fter a fueling event is set by the binding 
energy of the bulge regardless of the fueling mechanism (|Younger et al. 2008). In order to buildup a mass as large as 
M B ! ft X , multiple fueling episodes would likely be necessary, as is seen in cosmological simulations (e.g. JLi et al.ll2007t 
Di Matteo et al. 2008; Siiacki ct al. 2009). However, it is unlikely that Mg^ x could be increased significantly beyond 
the observed value as additional fueling mechanisms would likely result in only a small increase in the binding energy 
of the bulge. Therefore, additional fueling episodes would immediately result in the SMBH unbinding the accreting 
gas, preventing significant additional growth. Indeed, if our estimated BLQSO lifetime oUbl ~ 150 Myr is correct for 
SMBHs with Mbh ~ 1O 9 M , and if the SMBH is radiating at a typical Eddington ratio of L/LEdd ~ 0.1 during the 
BLQSO phase, then it would only accrete an additional ~ 10 8 M Q , assuming a radiative efficiency of e r = 0.1. If the 
distribution of Eddington ratios is a power-law with p(TEdd) oc T^^ 5 , as suggested by the discussion in § 15 . 1 . H then 
the typical Eddington ratio would be L/LEdd 0.1, suggesting even less growth during the BLQSO phase. Moreover, 
additional growth via black hole mergers, as might be expected from 'dry' mergers of galaxies, is also unlikely to lead 
to significant additional growth, as the mass of the additional SMBH is likely to be negligible compared to Mgg . And 
finally, additional growth t hrough radiatiyely inefficient low -accretion rate modes does not contribute significantly to 
the black hole's final mass (jHopkins et al.ll2006d: iCaol 120071 ). Therefore, our estimated value of Mgg* - 3 x 10 10 M Q 
should be representative of the most massive SMBH for both active and inactive SMBHs. 

Our results are consistent with previous data-based work that has attempted to map the distribution and growth 
of SMBHs. Most previous observational work has mapped SMBH growth by using the Mbh-& relationship to in- 
fer the local BHMF, and then stepped bac kward in time using the quasar luminosity function to infer the con- 
tribu t ion to the BHMF from accretion (e.g., iSoltanl H98l lYu fc Tremaindl2002t iShankar et "all 120041 : iMarconi et"all 
120041 : iMerloni fc Heinj I2008D . In ad dition, there have been attempts to predict the BHMF and its evolution for 
all S MBHs, or all a ctive SMBHs_(e.g.. lVolonteri et alj|2003t ISiiacki et al]|2007t iHopkins et aLll2008at iDi Matteo et al.l 
20081 iSiiacki et al.ll2009HShenll2009ft . While not directly comparable to these studies, as we focus on broad line quasars, 
our results and conclusions are qualitatively consistent with previous observational and theoretical work in that we 
find evidence for self-regulated SMBH growth, black hole downsizing, and BLQSO lifetimes of t B L ~ 150 Myr. 

6. SUMMARY 

Our main results are: 

• We have, for the first time, obtained an estimate of the black hole mass function for broad-line quasars that 
self-consistently corrects for incompleteness and the statistical uncertainty in the mass estimates derived from 
the broad emission lines in a statistically rigorous manner. Our estimated BHMF was obtained using data from 
the SDSS DR3 quasar sample. 

• The standard deviation in the statistical error of the broad line mass estimates is less than the commonly used 
~ 0.4 dex within the range of luminosity and redshift probed in our analysis. This may be due to correlation 
between the error in the mass estimates and luminosity and/or redshift, or a dependence of the standard deviation 
of the error on luminosity and/or redshift. When we treat the standard deviation in the statistical error as a 
free parameter, we estimate that the the Mg Il-based mass estimates scatter about the reverberation mapping 
estimates with an amplitude of » 0.18 dex, while the C IV-based estimate scatter with an amplitude of « 0.13 
dex. 
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• We find evidence for cosmic downsizing among BLQSOs, where the number density of BLQSOs peaks at higher 
redshift with increasing black hole mass. 

• We find that the comoving mass density of SMBHs in BLQSOs peaks at z ~ 2. We use our estimate for the 
BHMF to place constraints on the duty cycle, S(Mbh> z ), an d lifetime, t B L, for BLQSOs. The duty cycle at z = 1 
is constrained to be S > 0.01 at M BH ~ 10 9 M©, falling to 8 > 10~ 5 at M BH ~ 10 10 M Q . We estimate the lifetime 
of the BLQSO phase for SMBHs of M BH = 10 9 M Q at z = 1 to be t BL = 150 ± 15 Myr. However, we will have 
underestimated the BLQSO lifetime if there is a population of M BB = 10 9 M Q SMBHs that never experience a 
BLQSO phase along our line of sight, or if the local number density of these black holes is significantly larger than 
the z = 1 number density of these black holes. We argue that our estimated BLQSO lifetime, in combination 
with the estimated Eddington ratio distribution, suggests that most of a SMBH's growth occurs when it is not 
seen as a BLQSO and accreting at a higher rate, and that BLQSO activity represents a short phase that most 
SMBHs go through, consistent with self-regulated growth models. 

• We estimate that the most massive SMBH that could be seen as a BLQSO is M BH w 3 x lO lo A/ . This SMBH 
would most likely be seen as a BLQSO at z > 2. While largely in agreement with previous work, we have for 
the first time obtained statistically rigorous constraints on the value of Mg^ x and its redshift. 

• Assuming a constant bolometric correction of C1350 = 4.3 (Vester gaard fe Osmer|[2009f) . our inferred distribution 
of Eddington ratios peaks at L/L,Edd ~ 0.05 and has a dispersion of ~ 0.4 dex. Compared to previous work, 
our inferred Eddington ratio distribution is broader and shifted toward lower values of L/L B dd, showing that 
previous estimated distributions of L/ L B dd were significantly affected by incompleteness. We therefore provide 
evidence that most BLQSOs are not radiating at or near the Eddington limit, and that there is a large dispersion 
in Eddington ratio for BLQSOs. In addition, we also find that the number density of BLQSOs increases steeply 
toward lower values of L/ Lsdd, consistent with models where the BLQSO phase occurs when the fuel supply is 
dwindling or halted. 

• The evolution of the cosmological SMBH mass density for BLQ SOs tracks th e evolution in the cosmological 
accretion rate density of SMBHs estimated from variations of the iSoltanl {1982) argument (jMarconi et alJ 120041 : 

iMerloni "fc Heinz 2008; Shankar ct al. 2009). This result, in combination with our estimated BLQSO lifetime and 
Eddington ratio distribution, are qualitatively consistent with models of self-regulated S MBH growth, with the 
BLQSO phase occuring at t he end of the SMBHs obscured Eddington-limited growth (e.g.. I Hopkins fc Hernquisd 

1200a IHopkins et al.ll2006bh . 
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